<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<title>AMBHAS: /home/tomer/svn/ambhas/ambhas/richards.py Source File</title>

<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
<link href="navtree.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="resize.js"></script>
<script type="text/javascript" src="navtree.js"></script>
<script type="text/javascript">
  $(document).ready(initResizable);
</script>
<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
  $(document).ready(function() { searchBox.OnSelectItem(0); });
</script>

</head>
<body>
<div id="top"><!-- do not remove this div! -->


<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  
  
  <td style="padding-left: 0.5em;">
   <div id="projectname">AMBHAS
   
   </div>
   
  </td>
  
  
  
   
   <td>        <div id="MSearchBox" class="MSearchBoxInactive">
        <span class="left">
          <img id="MSearchSelect" src="search/mag_sel.png"
               onmouseover="return searchBox.OnSearchSelectShow()"
               onmouseout="return searchBox.OnSearchSelectHide()"
               alt=""/>
          <input type="text" id="MSearchField" value="Search" accesskey="S"
               onfocus="searchBox.OnSearchFieldFocus(true)" 
               onblur="searchBox.OnSearchFieldFocus(false)" 
               onkeyup="searchBox.OnSearchFieldChange(event)"/>
          </span><span class="right">
            <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
          </span>
        </div>
</td>
   
  
 </tr>
 </tbody>
</table>
</div>

<!-- Generated by Doxygen 1.7.6.1 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
</div>
<div id="side-nav" class="ui-resizable side-nav-resizable">
  <div id="nav-tree">
    <div id="nav-tree-contents">
    </div>
  </div>
  <div id="splitbar" style="-moz-user-select:none;" 
       class="ui-resizable-handle">
  </div>
</div>
<script type="text/javascript">
  initNavTree('richards_8py.html','');
</script>
<div id="doc-content">
<div class="header">
  <div class="headertitle">
<div class="title">richards.py</div>  </div>
</div><!--header-->
<div class="contents">
<a href="richards_8py.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a><a class="code" href="namespaceambhas_1_1richards.html">00001</a> <span class="comment"># -*- coding: utf-8 -*-</span>
<a name="l00002"></a>00002 <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00003"></a>00003 <span class="stringliteral">Created on Mon Mar 12 17:41:54 2012</span>
<a name="l00004"></a>00004 <span class="stringliteral"></span>
<a name="l00005"></a>00005 <span class="stringliteral">@author: sat kumar tomer</span>
<a name="l00006"></a>00006 <span class="stringliteral">@email: satkumartomer@gmail.com</span>
<a name="l00007"></a>00007 <span class="stringliteral">@website: www.ambhas.com</span>
<a name="l00008"></a>00008 <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00009"></a>00009 
<a name="l00010"></a>00010 <span class="keyword">from</span> __future__ <span class="keyword">import</span> division
<a name="l00011"></a>00011 <span class="keyword">import</span> numpy <span class="keyword">as</span> np
<a name="l00012"></a>00012 <span class="keyword">import</span> xlrd
<a name="l00013"></a>00013 <span class="keyword">from</span> Scientific.IO <span class="keyword">import</span> NetCDF <span class="keyword">as</span> nc
<a name="l00014"></a>00014 <span class="keyword">import</span> datetime
<a name="l00015"></a>00015 <span class="keyword">import</span> matplotlib.pyplot <span class="keyword">as</span> plt
<a name="l00016"></a>00016 <span class="keyword">from</span> BIP.Bayes.lhs <span class="keyword">import</span> lhs
<a name="l00017"></a>00017 <span class="keyword">from</span> scipy <span class="keyword">import</span> stats
<a name="l00018"></a>00018 <span class="comment">#np.seterr(all=&#39;raise&#39;)</span>
<a name="l00019"></a>00019 
<a name="l00020"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html">00020</a> <span class="keyword">class </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html">RICHARDS_1D</a>():
<a name="l00021"></a>00021     <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00022"></a>00022 <span class="stringliteral">    This is the main class of the RICHARDS_1D.</span>
<a name="l00023"></a>00023 <span class="stringliteral">    This simulates the flow in unsaturated porus media</span>
<a name="l00024"></a>00024 <span class="stringliteral">    </span>
<a name="l00025"></a>00025 <span class="stringliteral">    This will read the input data,</span>
<a name="l00026"></a>00026 <span class="stringliteral">    do the processing</span>
<a name="l00027"></a>00027 <span class="stringliteral">    and then write the output files</span>
<a name="l00028"></a>00028 <span class="stringliteral">    </span>
<a name="l00029"></a>00029 <span class="stringliteral">    &quot;&quot;&quot;</span>
<a name="l00030"></a>00030      
<a name="l00031"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#abe0e4d7b40c9db697f847bd46129d7d0">00031</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#abe0e4d7b40c9db697f847bd46129d7d0">__init__</a>(self,input_file):
<a name="l00032"></a>00032         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00033"></a>00033 <span class="stringliteral">        Input:</span>
<a name="l00034"></a>00034 <span class="stringliteral">            input_file: the file which contains all the information</span>
<a name="l00035"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">00035</a> <span class="stringliteral">            including forcing and parameters.</span>
<a name="l00036"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f">00036</a> <span class="stringliteral">        &quot;&quot;&quot;</span>      
<a name="l00037"></a>00037         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a> = input_file
<a name="l00038"></a>00038         
<a name="l00039"></a>00039         <span class="comment"># read the input data</span>
<a name="l00040"></a>00040         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a543acf445325c80ee0602aa414f472c1">_read_input</a>()
<a name="l00041"></a>00041         
<a name="l00042"></a>00042         <span class="comment"># initialize the variables and output file</span>
<a name="l00043"></a>00043         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a318ed5c6afbde685776cacdb4284c068">initialize</a>()
<a name="l00044"></a>00044         
<a name="l00045"></a>00045         <span class="comment">################ run the model ########################</span>
<a name="l00046"></a>00046         <span class="keywordflow">for</span> t <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2702f64024162cdafa67eddfe7dc016d">max_t</a>):
<a name="l00047"></a>00047             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a> = t
<a name="l00048"></a>00048               
<a name="l00049"></a>00049             <span class="comment"># get forcing data at current time step        </span>
<a name="l00050"></a>00050             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a5b9ee80409ee56964194b5892ded099d">_get_forcing</a>()
<a name="l00051"></a>00051             
<a name="l00052"></a>00052             <span class="comment"># call the unsat module</span>
<a name="l00053"></a>00053             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ad9340e6ce76db76ed3dc7ed6e855584f">_unsat</a>()
<a name="l00054"></a>00054                 
<a name="l00055"></a>00055         self.nc_file.close() <span class="comment"># close the output file</span>
<a name="l00056"></a>00056          
<a name="l00057"></a>00057      
<a name="l00058"></a>00058             
<a name="l00059"></a>00059     <span class="keyword">def </span>_read_input(self):
<a name="l00060"></a>00060         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00061"></a>00061 <span class="stringliteral">        This checks if all the required input sheets are present in the xls file,</span>
<a name="l00062"></a>00062 <span class="stringliteral">        read the data from input file, which can be used later in other functions</span>
<a name="l00063"></a>00063 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00064"></a>00064     
<a name="l00065"></a>00065         <span class="comment"># list of required files in the input directory</span>
<a name="l00066"></a>00066         input_sheets = [<span class="stringliteral">&#39;ind&#39;</span>, <span class="stringliteral">&#39;forcing&#39;</span>, <span class="stringliteral">&#39;initial_condition&#39;</span>, <span class="stringliteral">&#39;units&#39;</span>, <span class="stringliteral">&#39;temporal_info&#39;</span>,
<a name="l00067"></a>00067                        <span class="stringliteral">&#39;spatial_info&#39;</span>, <span class="stringliteral">&#39;soil_hyd_par&#39;</span>, <span class="stringliteral">&#39;output_par&#39;</span>]
<a name="l00068"></a>00068         
<a name="l00069"></a>00069         <span class="comment"># check if all the required sheets are present or not</span>
<a name="l00070"></a>00070         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aa38dc52915838e2e166b1ee1af8a34ef">_check_sheets</a>(input_sheets, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00071"></a>00071         
<a name="l00072"></a>00072         <span class="comment"># read the legend</span>
<a name="l00073"></a>00073         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a3ecc113a0ee750b7700cfca9215e5a97">_read_ind</a>()
<a name="l00074"></a>00074         
<a name="l00075"></a>00075         <span class="comment"># read the spatial data</span>
<a name="l00076"></a>00076         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a3fb57f5fc2b5bc190da51b79a9677755">_read_spatial</a>()
<a name="l00077"></a>00077         
<a name="l00078"></a>00078         <span class="comment"># read the temporal data</span>
<a name="l00079"></a>00079         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a7ea1b0013ddb45271fa721a677be1568">_read_temporal</a>()
<a name="l00080"></a>00080 
<a name="l00081"></a>00081         <span class="comment"># read the units </span>
<a name="l00082"></a>00082         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a01b38370e7172080e015734328b6f6c2">_read_units</a>()
<a name="l00083"></a>00083         
<a name="l00084"></a>00084         <span class="comment"># read the initial condition</span>
<a name="l00085"></a>00085         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aa7d5b17074db9d3f27dd9262581c892d">_read_initial_condition</a>()
<a name="l00086"></a>00086         
<a name="l00087"></a>00087         <span class="comment"># read the soil hydraulic properties data</span>
<a name="l00088"></a>00088         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a8f90854a5f37f92a3a4bf5c30b9201b7">_read_shp</a>()
<a name="l00089"></a>00089         
<a name="l00090"></a>00090         <span class="comment"># read the forcing infomation</span>
<a name="l00091"></a>00091         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a1433858e28ba363994050df463f6c0d0">_read_forcing</a>()
<a name="l00092"></a>00092         
<a name="l00093"></a>00093         <span class="comment"># read the outfile name</span>
<a name="l00094"></a>00094         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a51f0263bd2da91bf3410ab1784519334">_read_ofile_name</a>()
<a name="l00095"></a>00095         
<a name="l00096"></a>00096         <span class="comment"># print the reading status</span>
<a name="l00097"></a>00097         output_message = <span class="stringliteral">&#39;Input data reading completed sucessfully&#39;</span>
<a name="l00098"></a>00098         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ab183d7511d72f3612a15ef2e20bdfc37">_colored_output</a>(output_message, 32)
<a name="l00099"></a>00099     
<a name="l00100"></a>00100     <span class="keyword">def </span>_check_sheets(self, check_sheets, check_file):
<a name="l00101"></a>00101         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00102"></a>00102 <span class="stringliteral">        This functions check if all the sheets needed to model are present  </span>
<a name="l00103"></a>00103 <span class="stringliteral">        in check_file</span>
<a name="l00104"></a>00104 <span class="stringliteral">        </span>
<a name="l00105"></a>00105 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00106"></a>00106         <span class="comment"># open the xls file and get its sheets</span>
<a name="l00107"></a>00107         foo = xlrd.open_workbook(check_file)
<a name="l00108"></a>00108         check_sheet_names = foo.sheet_names()
<a name="l00109"></a>00109         
<a name="l00110"></a>00110         <span class="keywordflow">for</span> check_sheet <span class="keywordflow">in</span> check_sheets:
<a name="l00111"></a>00111             <span class="keywordflow">if</span> check_sheet <span class="keywordflow">not</span> <span class="keywordflow">in</span> check_sheet_names:
<a name="l00112"></a>00112                 output_message = <span class="stringliteral">&#39;The sheet &#39;</span> + check_sheet + <span class="stringliteral">&#39; is missing&#39;</span>
<a name="l00113"></a>00113                 self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ab183d7511d72f3612a15ef2e20bdfc37">_colored_output</a>(output_message,31)
<a name="l00114"></a>00114             
<a name="l00115"></a>00115 
<a name="l00116"></a>00116     <span class="keyword">def </span>_read_ind(self):
<a name="l00117"></a>00117         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00118"></a>00118 <span class="stringliteral">        Read the ind sheet</span>
<a name="l00119"></a>00119 <span class="stringliteral">        legend stores the information about the indices of other properties,</span>
<a name="l00120"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">00120</a> <span class="stringliteral">        which would be used by all other properties reading functions</span>
<a name="l00121"></a>00121 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00122"></a>00122         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00123"></a>00123         sheet = book.sheet_by_name(<span class="stringliteral">&#39;ind&#39;</span>)
<a name="l00124"></a>00124         <span class="comment"># dont read the first line of the xls file</span>
<a name="l00125"></a>00125         ind = {}
<a name="l00126"></a>00126         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(sheet.nrows-1):
<a name="l00127"></a>00127             ind[str(sheet.cell_value(i+1,0))] = int(sheet.cell_value(i+1,1))
<a name="l00128"></a>00128                 
<a name="l00129"></a>00129         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">ind</a> = ind
<a name="l00130"></a>00130 
<a name="l00131"></a>00131     <span class="keyword">def </span>_read_spatial(self):
<a name="l00132"></a>00132         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00133"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a98e65f23eb8575aa56b137588c4f5a50">00133</a> <span class="stringliteral">        Read the spatial info</span>
<a name="l00134"></a>00134 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00135"></a>00135         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00136"></a>00136         sheet = book.sheet_by_name(<span class="stringliteral">&#39;spatial_info&#39;</span>)
<a name="l00137"></a>00137         <span class="comment"># get the row number from the ind</span>
<a name="l00138"></a>00138         j = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">ind</a>[<span class="stringliteral">&#39;spatial_info&#39;</span>]
<a name="l00139"></a>00139         no_layer = int(sheet.cell_value(j,1))
<a name="l00140"></a>00140         dz = sheet.cell_value(j,2)
<a name="l00141"></a>00141                 
<a name="l00142"></a>00142         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a98e65f23eb8575aa56b137588c4f5a50">no_layer</a> = no_layer
<a name="l00143"></a>00143         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a54a1bf0c73e3c943a5aadf73751e058b">dz</a> = dz
<a name="l00144"></a>00144     
<a name="l00145"></a>00145     <span class="keyword">def </span>_read_temporal(self):
<a name="l00146"></a>00146         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00147"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a4883d42564004396dda7e720bdbeda86">00147</a> <span class="stringliteral">        Read the temporal info</span>
<a name="l00148"></a>00148 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00149"></a>00149         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00150"></a>00150         sheet = book.sheet_by_name(<span class="stringliteral">&#39;temporal_info&#39;</span>)
<a name="l00151"></a>00151         <span class="comment">#get the row number from the ind</span>
<a name="l00152"></a>00152         j = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">ind</a>[<span class="stringliteral">&#39;temporal_info&#39;</span>]
<a name="l00153"></a>00153         dt = sheet.cell_value(j,1)
<a name="l00154"></a>00154         final_time = sheet.cell_value(j,2)
<a name="l00155"></a>00155         
<a name="l00156"></a>00156         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57c5d72846104a31fa67815aa57f700d">dt_flux</a> = dt
<a name="l00157"></a>00157         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a4883d42564004396dda7e720bdbeda86">final_time</a> = final_time
<a name="l00158"></a>00158     
<a name="l00159"></a>00159     
<a name="l00160"></a>00160     <span class="keyword">def </span>_read_units(self):
<a name="l00161"></a>00161         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00162"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#afd9d919140433fea93c16444765e2b20">00162</a> <span class="stringliteral">        read the units of the forcing data</span>
<a name="l00163"></a>00163 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00164"></a>00164         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00165"></a>00165         sheet = book.sheet_by_name(<span class="stringliteral">&#39;units&#39;</span>)
<a name="l00166"></a>00166         <span class="comment">#get the row number from the ind</span>
<a name="l00167"></a>00167         j = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">ind</a>[<span class="stringliteral">&#39;units&#39;</span>]
<a name="l00168"></a>00168         forcing_units = {}
<a name="l00169"></a>00169         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(sheet.ncols-1):
<a name="l00170"></a>00170             forcing_units[str(sheet.cell_value(0,i+1))] = str(sheet.cell_value(j,i+1))
<a name="l00171"></a>00171         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#afd9d919140433fea93c16444765e2b20">forcing_units</a> = forcing_units
<a name="l00172"></a>00172     
<a name="l00173"></a>00173     <span class="keyword">def </span>_read_initial_condition(self):
<a name="l00174"></a>00174         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00175"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ade9c77173908940774a8581b18cd7a74">00175</a> <span class="stringliteral">        read initial condition</span>
<a name="l00176"></a>00176 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00177"></a>00177         <span class="comment">#get the row number from the ind</span>
<a name="l00178"></a>00178         j = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">ind</a>[<span class="stringliteral">&#39;initial_condition&#39;</span>]
<a name="l00179"></a>00179         
<a name="l00180"></a>00180         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00181"></a>00181         sheet = book.sheet_by_name(<span class="stringliteral">&#39;initial_condition&#39;</span>)
<a name="l00182"></a>00182         theta_0 = sheet.cell_value(j,1)
<a name="l00183"></a>00183         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ade9c77173908940774a8581b18cd7a74">theta</a> = np.tile(theta_0,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a98e65f23eb8575aa56b137588c4f5a50">no_layer</a>)
<a name="l00184"></a>00184         
<a name="l00185"></a>00185     <span class="keyword">def </span>_read_shp(self):
<a name="l00186"></a>00186         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00187"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">00187</a> <span class="stringliteral">        read the soil hydraulic parameters</span>
<a name="l00188"></a>00188 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00189"></a>00189         <span class="comment">#get the row number from the ind</span>
<a name="l00190"></a>00190         j = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">ind</a>[<span class="stringliteral">&#39;soil_hyd_par&#39;</span>]
<a name="l00191"></a>00191         
<a name="l00192"></a>00192         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00193"></a>00193         sheet = book.sheet_by_name(<span class="stringliteral">&#39;soil_hyd_par&#39;</span>)
<a name="l00194"></a>00194         soil_par = {}
<a name="l00195"></a>00195         soil_par[<span class="stringliteral">&#39;thetar&#39;</span>] = sheet.cell_value(j,1)
<a name="l00196"></a>00196         soil_par[<span class="stringliteral">&#39;thetas&#39;</span>] = sheet.cell_value(j,2)
<a name="l00197"></a>00197         soil_par[<span class="stringliteral">&#39;alpha&#39;</span>] = sheet.cell_value(j,3)
<a name="l00198"></a>00198         soil_par[<span class="stringliteral">&#39;n&#39;</span>] = sheet.cell_value(j,4)
<a name="l00199"></a>00199         soil_par[<span class="stringliteral">&#39;Ks&#39;</span>] = sheet.cell_value(j,5)
<a name="l00200"></a>00200         soil_par[<span class="stringliteral">&#39;l&#39;</span>] = sheet.cell_value(j,6)
<a name="l00201"></a>00201         soil_par[<span class="stringliteral">&#39;evap_0&#39;</span>] = sheet.cell_value(j,7)
<a name="l00202"></a>00202         soil_par[<span class="stringliteral">&#39;evap_1&#39;</span>] = sheet.cell_value(j,8)
<a name="l00203"></a>00203         soil_par[<span class="stringliteral">&#39;m&#39;</span>] = 1-1/soil_par[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l00204"></a>00204         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a> = soil_par
<a name="l00205"></a>00205     
<a name="l00206"></a>00206     
<a name="l00207"></a>00207     <span class="keyword">def </span>_read_forcing(self):
<a name="l00208"></a>00208         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00209"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a071af44b049e83375c32281d32e2649e">00209</a> <span class="stringliteral">        read the forcing data from xls file</span>
<a name="l00210"></a>00210 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00211"></a>00211         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00212"></a>00212         sheet = book.sheet_by_name(<span class="stringliteral">&#39;forcing&#39;</span>)
<a name="l00213"></a>00213         
<a name="l00214"></a>00214         data_len = sheet.nrows-1
<a name="l00215"></a>00215         year = np.zeros(data_len)
<a name="l00216"></a>00216         doy = np.zeros(data_len)
<a name="l00217"></a>00217         rain = np.zeros(data_len)
<a name="l00218"></a>00218         pet = np.zeros(data_len)
<a name="l00219"></a>00219             
<a name="l00220"></a>00220         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> xrange(data_len):
<a name="l00221"></a>00221             year[i] = sheet.cell_value(i+1,0)
<a name="l00222"></a>00222             doy[i] = sheet.cell_value(i+1,1)
<a name="l00223"></a>00223             rain[i] = sheet.cell_value(i+1,2)
<a name="l00224"></a>00224             pet[i] = sheet.cell_value(i+1,3)
<a name="l00225"></a>00225                     
<a name="l00226"></a>00226         
<a name="l00227"></a>00227         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a071af44b049e83375c32281d32e2649e">year</a> = year
<a name="l00228"></a>00228         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a928a3af97f3ed434e2410cf5ea989f86">doy</a> = doy
<a name="l00229"></a>00229         
<a name="l00230"></a>00230         <span class="comment"># if forcing data was in mm units, covert into m</span>
<a name="l00231"></a>00231         <span class="keywordflow">if</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#afd9d919140433fea93c16444765e2b20">forcing_units</a>[<span class="stringliteral">&#39;rain&#39;</span>] == <span class="stringliteral">&#39;mm&#39;</span>:
<a name="l00232"></a>00232             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a6201ea018543affc7ddfefe044e9c581">rain</a> = rain/1000.0
<a name="l00233"></a>00233         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#afd9d919140433fea93c16444765e2b20">forcing_units</a>[<span class="stringliteral">&#39;rain&#39;</span>] == <span class="stringliteral">&#39;m&#39;</span>:
<a name="l00234"></a>00234             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a6201ea018543affc7ddfefe044e9c581">rain</a> = rain
<a name="l00235"></a>00235         <span class="keywordflow">else</span>:
<a name="l00236"></a>00236             <span class="keywordflow">raise</span> ValueError(<span class="stringliteral">&quot;The units of rain should be either &#39;mm&#39; or &#39;m&#39; &quot;</span>)
<a name="l00237"></a>00237 
<a name="l00238"></a>00238         <span class="keywordflow">if</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#afd9d919140433fea93c16444765e2b20">forcing_units</a>[<span class="stringliteral">&#39;pet&#39;</span>] == <span class="stringliteral">&#39;mm&#39;</span>:
<a name="l00239"></a>00239             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2d0fb0fdaba97b19ea0d2ce0cb06d443">pet</a> = pet/1000.0
<a name="l00240"></a>00240         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#afd9d919140433fea93c16444765e2b20">forcing_units</a>[<span class="stringliteral">&#39;pet&#39;</span>] == <span class="stringliteral">&#39;m&#39;</span>:
<a name="l00241"></a>00241             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2d0fb0fdaba97b19ea0d2ce0cb06d443">pet</a> = pet
<a name="l00242"></a>00242         <span class="keywordflow">else</span>:
<a name="l00243"></a>00243             <span class="keywordflow">raise</span> ValueError(<span class="stringliteral">&quot;The units of PET should be either &#39;mm&#39; or &#39;m&#39; &quot;</span>)
<a name="l00244"></a>00244             
<a name="l00245"></a>00245       
<a name="l00246"></a>00246 
<a name="l00247"></a>00247     <span class="keyword">def </span>_read_ofile_name(self):
<a name="l00248"></a>00248         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00249"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2bd6df8896bc909494971fe14bed36ea">00249</a> <span class="stringliteral">        read the forcing data from xls file</span>
<a name="l00250"></a>00250 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00251"></a>00251         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a9116342a7ec04c2f655bd6991acd5a87">input_file</a>)
<a name="l00252"></a>00252         sheet = book.sheet_by_name(<span class="stringliteral">&#39;output_par&#39;</span>)
<a name="l00253"></a>00253         j = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef9169c10725dd24c126618b822b5442">ind</a>[<span class="stringliteral">&#39;output_par&#39;</span>]
<a name="l00254"></a>00254         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2bd6df8896bc909494971fe14bed36ea">ofile_name</a> = str(sheet.cell_value(j,1))
<a name="l00255"></a>00255 
<a name="l00256"></a>00256 
<a name="l00257"></a>00257     <span class="keyword">def </span>_colored_output(self, output_message, color):
<a name="l00258"></a>00258         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00259"></a>00259 <span class="stringliteral">        This functions print  the output_message in the color</span>
<a name="l00260"></a>00260 <span class="stringliteral">        Input:</span>
<a name="l00261"></a>00261 <span class="stringliteral">            output_messgae: the text you want to print</span>
<a name="l00262"></a>00262 <span class="stringliteral">            color: the color in which you want to print text, it could be one of:</span>
<a name="l00263"></a>00263 <span class="stringliteral">                30: Gray</span>
<a name="l00264"></a>00264 <span class="stringliteral">                31: Red</span>
<a name="l00265"></a>00265 <span class="stringliteral">                32: Green</span>
<a name="l00266"></a>00266 <span class="stringliteral">                33: Yellow</span>
<a name="l00267"></a>00267 <span class="stringliteral">                34: Blue</span>
<a name="l00268"></a>00268 <span class="stringliteral">                35: Magneta</span>
<a name="l00269"></a>00269 <span class="stringliteral">                66: Cyan</span>
<a name="l00270"></a>00270 <span class="stringliteral">                37: White</span>
<a name="l00271"></a>00271 <span class="stringliteral">                38: Crimson</span>
<a name="l00272"></a>00272 <span class="stringliteral">                41: Highlighted Red</span>
<a name="l00273"></a>00273 <span class="stringliteral">                42: Highlighted Green </span>
<a name="l00274"></a>00274 <span class="stringliteral">                43: Highlighted Brown </span>
<a name="l00275"></a>00275 <span class="stringliteral">                44: Highlighted Blue </span>
<a name="l00276"></a>00276 <span class="stringliteral">                45: Highlighted Magenta </span>
<a name="l00277"></a>00277 <span class="stringliteral">                46: Highlighted Cyan</span>
<a name="l00278"></a>00278 <span class="stringliteral">                47: Highlighted Gray </span>
<a name="l00279"></a>00279 <span class="stringliteral">                48: Highlighted Crimson </span>
<a name="l00280"></a>00280 <span class="stringliteral">        Output:</span>
<a name="l00281"></a>00281 <span class="stringliteral">            This returns None, but print the output in python shell</span>
<a name="l00282"></a>00282 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00283"></a>00283                 
<a name="l00284"></a>00284         print((<span class="stringliteral">&quot;\033[31m&quot;</span> +output_message+ <span class="stringliteral">&quot;\033[0m&quot;</span>).replace(<span class="stringliteral">&#39;31&#39;</span>,str(color)))
<a name="l00285"></a>00285 
<a name="l00286"></a>00286     <span class="keyword">def </span>_get_forcing(self):
<a name="l00287"></a>00287         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00288"></a>00288 <span class="stringliteral">        this will give the forcing at time t</span>
<a name="l00289"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a05f361eeb6230c24a1d7718395c71b3b">00289</a> <span class="stringliteral">        forcing are given in terms of L/T</span>
<a name="l00290"></a>00290 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00291"></a>00291         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a05f361eeb6230c24a1d7718395c71b3b">rain_cur</a> = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a6201ea018543affc7ddfefe044e9c581">rain</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>]/self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57c5d72846104a31fa67815aa57f700d">dt_flux</a>
<a name="l00292"></a>00292         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a6e3b8c0a4f36c24d8af1238d81c33945">pet_cur</a> = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2d0fb0fdaba97b19ea0d2ce0cb06d443">pet</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>]/self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57c5d72846104a31fa67815aa57f700d">dt_flux</a>
<a name="l00293"></a>00293                 
<a name="l00294"></a>00294         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a4473319b3c1ac79c3a1f65bb59334506">cur_year</a> = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a071af44b049e83375c32281d32e2649e">year</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>]
<a name="l00295"></a>00295         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a5dd83c846f2ef3fa7d7a3656e85df861">cur_doy</a> = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a928a3af97f3ed434e2410cf5ea989f86">doy</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>]
<a name="l00296"></a>00296         
<a name="l00297"></a>00297     
<a name="l00298"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aa83c63285e5f5efad117686fc9683af5">00298</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aa83c63285e5f5efad117686fc9683af5">smcf</a>(self, theta, thetar, thetas, alpha, m, n):
<a name="l00299"></a>00299         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00300"></a>00300 <span class="stringliteral">        smcf: calculate the smc</span>
<a name="l00301"></a>00301 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00302"></a>00302         Se = (theta-thetar)/(thetas-thetar)
<a name="l00303"></a>00303         Se[Se&lt;=0] = 0.0001
<a name="l00304"></a>00304         Se[Se&gt;=1] = 0.9999
<a name="l00305"></a>00305         <span class="keywordflow">try</span>:
<a name="l00306"></a>00306             smc = alpha*(thetas-thetar)*m*n*pow(Se,1/m+1)*pow(pow(Se,-1/m)-1,m)
<a name="l00307"></a>00307         <span class="keywordflow">except</span>:
<a name="l00308"></a>00308             <span class="keywordflow">print</span> thetas,thetar,m,n,Se,
<a name="l00309"></a>00309         <span class="keywordflow">return</span> smc
<a name="l00310"></a>00310     
<a name="l00311"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a0445a3a4323788bfeb39b0dd4bf22e4a">00311</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a0445a3a4323788bfeb39b0dd4bf22e4a">theta2psi</a>(self,theta, thetar, thetas, m, n, alpha):
<a name="l00312"></a>00312         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00313"></a>00313 <span class="stringliteral">        theta2psi: given the theta calculate the psi</span>
<a name="l00314"></a>00314 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00315"></a>00315         Se = (theta-thetar)/(thetas-thetar)
<a name="l00316"></a>00316         Se[Se&lt;=0] = 0.0001
<a name="l00317"></a>00317         Se[Se&gt;=1] = 0.9999
<a name="l00318"></a>00318         psi = -(1/alpha)*pow(pow(Se,-1/m)-1,1/n)
<a name="l00319"></a>00319         psi[psi&lt;-1e6] = -1e6
<a name="l00320"></a>00320         <span class="keywordflow">return</span> psi
<a name="l00321"></a>00321         
<a name="l00322"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef469b323d85eeaaa6ca6f1bcc5de350">00322</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aef469b323d85eeaaa6ca6f1bcc5de350">psi2theta</a>(self,psi, thetar, thetas, alpha, m, n):
<a name="l00323"></a>00323          <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00324"></a>00324 <span class="stringliteral">         psi2theta: given the theta calculate the pressure head</span>
<a name="l00325"></a>00325 <span class="stringliteral">         &quot;&quot;&quot;</span>
<a name="l00326"></a>00326          <span class="keywordflow">if</span> (psi&gt;=0):
<a name="l00327"></a>00327              theta = thetas
<a name="l00328"></a>00328          <span class="keywordflow">elif</span> psi&lt;-1e6:
<a name="l00329"></a>00329              theta = 1.01*thetar
<a name="l00330"></a>00330          <span class="keywordflow">else</span>:
<a name="l00331"></a>00331              theta = thetar+(thetas-thetar)*pow(1+pow(abs(alpha*psi),n),-m)
<a name="l00332"></a>00332          <span class="keywordflow">return</span> theta
<a name="l00333"></a>00333          
<a name="l00334"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57edd92079beb2139337a15db7868836">00334</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57edd92079beb2139337a15db7868836">theta2kr</a>(self,theta, thetar, thetas, m, l, Ks):
<a name="l00335"></a>00335         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00336"></a>00336 <span class="stringliteral">        theta2kr: given the theta, calculate the kr </span>
<a name="l00337"></a>00337 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00338"></a>00338         Se = (theta-thetar)/(thetas-thetar)
<a name="l00339"></a>00339         Se[Se&lt;0] = 0.00001
<a name="l00340"></a>00340         Se[Se&gt;1] = 0.99999
<a name="l00341"></a>00341         kr = Ks*(pow(Se,l))*pow(1-pow(1-pow(Se,1/m),m),2)
<a name="l00342"></a>00342         kr[Se&lt;0] = 0
<a name="l00343"></a>00343         kr[Se&gt;1] = Ks
<a name="l00344"></a>00344         
<a name="l00345"></a>00345         <span class="keywordflow">return</span> kr
<a name="l00346"></a>00346     
<a name="l00347"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a318ed5c6afbde685776cacdb4284c068">00347</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a318ed5c6afbde685776cacdb4284c068">initialize</a>(self):
<a name="l00348"></a>00348         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00349"></a>00349 <span class="stringliteral">        this initializes all the required variables</span>
<a name="l00350"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a60c00cb4c4e4c0934eec2b046b970e72">00350</a> <span class="stringliteral">        and open the netcdf file for writting</span>
<a name="l00351"></a>00351 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00352"></a>00352         max_t = int(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a4883d42564004396dda7e720bdbeda86">final_time</a>/self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57c5d72846104a31fa67815aa57f700d">dt_flux</a>)
<a name="l00353"></a>00353         <span class="comment">#max_t = 56</span>
<a name="l00354"></a>00354         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2702f64024162cdafa67eddfe7dc016d">max_t</a> = max_t
<a name="l00355"></a>00355         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ae87ba3c49ae83bcfba2b17b88bcff452">iter_dt</a> = 1
<a name="l00356"></a>00356                         
<a name="l00357"></a>00357         <span class="comment"># open file for writing</span>
<a name="l00358"></a>00358         file = nc.NetCDFFile(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2bd6df8896bc909494971fe14bed36ea">ofile_name</a>, <span class="stringliteral">&#39;w&#39;</span>)
<a name="l00359"></a>00359         setattr(file, <span class="stringliteral">&#39;title&#39;</span>, <span class="stringliteral">&#39;output of the model ambhas.richards&#39;</span>)
<a name="l00360"></a>00360         now = datetime.datetime.now()
<a name="l00361"></a>00361         setattr(file, <span class="stringliteral">&#39;description&#39;</span>, <span class="stringliteral">&#39;The model was run at %s&#39;</span>%(now.ctime()))
<a name="l00362"></a>00362         file.createDimension(<span class="stringliteral">&#39;depth&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a98e65f23eb8575aa56b137588c4f5a50">no_layer</a>)
<a name="l00363"></a>00363         file.createDimension(<span class="stringliteral">&#39;time&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2702f64024162cdafa67eddfe7dc016d">max_t</a>+1)
<a name="l00364"></a>00364         
<a name="l00365"></a>00365         <span class="comment"># depth</span>
<a name="l00366"></a>00366         varDims = <span class="stringliteral">&#39;depth&#39;</span>,
<a name="l00367"></a>00367         depth = file.createVariable(<span class="stringliteral">&#39;depth&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00368"></a>00368         depth.units = <span class="stringliteral">&#39;m&#39;</span>
<a name="l00369"></a>00369         depth[:] = np.tile(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a54a1bf0c73e3c943a5aadf73751e058b">dz</a>,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a98e65f23eb8575aa56b137588c4f5a50">no_layer</a>).cumsum()-self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a54a1bf0c73e3c943a5aadf73751e058b">dz</a>/2
<a name="l00370"></a>00370         
<a name="l00371"></a>00371         <span class="comment"># time (year and doy)</span>
<a name="l00372"></a>00372         varDims = <span class="stringliteral">&#39;time&#39;</span>,
<a name="l00373"></a>00373         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a60c00cb4c4e4c0934eec2b046b970e72">nc_year</a> = file.createVariable(<span class="stringliteral">&#39;year&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00374"></a>00374         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a8040d34cfb0fec5d81c37d8be3e0c569">nc_doy</a> = file.createVariable(<span class="stringliteral">&#39;doy&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00375"></a>00375         
<a name="l00376"></a>00376         <span class="comment"># soil moisture</span>
<a name="l00377"></a>00377         varDims = <span class="stringliteral">&#39;depth&#39;</span>,<span class="stringliteral">&#39;time&#39;</span>
<a name="l00378"></a>00378         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#af7a9b280883bbb1ced757befda5155a3">nc_sm</a> = file.createVariable(<span class="stringliteral">&#39;sm&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00379"></a>00379         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#af7a9b280883bbb1ced757befda5155a3">nc_sm</a>.units = <span class="stringliteral">&#39;v/v&#39;</span>
<a name="l00380"></a>00380         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#af7a9b280883bbb1ced757befda5155a3">nc_sm</a>[:,0] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ade9c77173908940774a8581b18cd7a74">theta</a>
<a name="l00381"></a>00381         
<a name="l00382"></a>00382         <span class="comment"># recharge and aet</span>
<a name="l00383"></a>00383         varDims = <span class="stringliteral">&#39;time&#39;</span>,
<a name="l00384"></a>00384         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#acdbad2c6ebdb951bf36b1177eb2d9023">nc_aet</a> = file.createVariable(<span class="stringliteral">&#39;aet&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00385"></a>00385         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#acdbad2c6ebdb951bf36b1177eb2d9023">nc_aet</a>.units = <span class="stringliteral">&#39;mm&#39;</span>
<a name="l00386"></a>00386         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#af46617063a50d73918d2f3b82163612e">nc_recharge</a> = file.createVariable(<span class="stringliteral">&#39;recharge&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00387"></a>00387         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#af46617063a50d73918d2f3b82163612e">nc_recharge</a>.units = <span class="stringliteral">&#39;mm&#39;</span>
<a name="l00388"></a>00388         
<a name="l00389"></a>00389         <span class="comment"># soil_par</span>
<a name="l00390"></a>00390         setattr(file, <span class="stringliteral">&#39;thetar&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;thetar&#39;</span>])
<a name="l00391"></a>00391         setattr(file, <span class="stringliteral">&#39;thetas&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;thetas&#39;</span>])
<a name="l00392"></a>00392         setattr(file, <span class="stringliteral">&#39;alpha&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;alpha&#39;</span>])
<a name="l00393"></a>00393         setattr(file, <span class="stringliteral">&#39;n&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;n&#39;</span>])
<a name="l00394"></a>00394         setattr(file, <span class="stringliteral">&#39;Ks&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;Ks&#39;</span>])
<a name="l00395"></a>00395         setattr(file, <span class="stringliteral">&#39;l&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;l&#39;</span>])        
<a name="l00396"></a>00396         
<a name="l00397"></a>00397         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a0b202186e73f83d8d81a27394dedeeef">nc_file</a> = file
<a name="l00398"></a>00398                 
<a name="l00399"></a>00399         
<a name="l00400"></a>00400     <span class="keyword">def </span>_unsat(self):
<a name="l00401"></a>00401         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00402"></a>00402 <span class="stringliteral">        top boundary: atmoshpheric</span>
<a name="l00403"></a>00403 <span class="stringliteral">        bottom boundary: gravity drainage</span>
<a name="l00404"></a>00404 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00405"></a>00405                
<a name="l00406"></a>00406         thetar = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;thetar&#39;</span>]
<a name="l00407"></a>00407         thetas = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;thetas&#39;</span>]
<a name="l00408"></a>00408         alpha = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;alpha&#39;</span>]
<a name="l00409"></a>00409         n = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l00410"></a>00410         m = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;m&#39;</span>]
<a name="l00411"></a>00411         l = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;l&#39;</span>]
<a name="l00412"></a>00412         Ks = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;Ks&#39;</span>]
<a name="l00413"></a>00413         nz = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a98e65f23eb8575aa56b137588c4f5a50">no_layer</a>
<a name="l00414"></a>00414                 
<a name="l00415"></a>00415         theta = 1.0*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ade9c77173908940774a8581b18cd7a74">theta</a>
<a name="l00416"></a>00416         
<a name="l00417"></a>00417         <span class="comment">#delta_theta = (np.abs(flux)).max()</span>
<a name="l00418"></a>00418         
<a name="l00419"></a>00419         iter_dt = max(24,int(np.ceil(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a05f361eeb6230c24a1d7718395c71b3b">rain_cur</a>*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57c5d72846104a31fa67815aa57f700d">dt_flux</a>*1000/0.15)))
<a name="l00420"></a>00420         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ae87ba3c49ae83bcfba2b17b88bcff452">iter_dt</a> = int(max(iter_dt,0.75*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ae87ba3c49ae83bcfba2b17b88bcff452">iter_dt</a>))
<a name="l00421"></a>00421         
<a name="l00422"></a>00422         <span class="comment">#if self.t == 56:</span>
<a name="l00423"></a>00423         <span class="comment">#    self.iter_dt = int(self.iter_dt*6)</span>
<a name="l00424"></a>00424         <span class="comment">#print self.iter_dt</span>
<a name="l00425"></a>00425         
<a name="l00426"></a>00426         recharge_day = 0
<a name="l00427"></a>00427         aet_day = 0
<a name="l00428"></a>00428         
<a name="l00429"></a>00429         <span class="comment"># check for time step</span>
<a name="l00430"></a>00430         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ae87ba3c49ae83bcfba2b17b88bcff452">iter_dt</a>):
<a name="l00431"></a>00431             dt = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57c5d72846104a31fa67815aa57f700d">dt_flux</a>/self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ae87ba3c49ae83bcfba2b17b88bcff452">iter_dt</a>
<a name="l00432"></a>00432             <span class="comment">#print dt</span>
<a name="l00433"></a>00433             <span class="comment"># top boundary value</span>
<a name="l00434"></a>00434             smi = (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ade9c77173908940774a8581b18cd7a74">theta</a>[0]-self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;evap_0&#39;</span>])/(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;evap_1&#39;</span>]-self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a556ff101921f19f742580bfe1ea8998e">soil_par</a>[<span class="stringliteral">&#39;evap_0&#39;</span>])
<a name="l00435"></a>00435             <span class="keywordflow">if</span> smi&lt;0: smi=0
<a name="l00436"></a>00436             <span class="keywordflow">if</span> smi&gt;1: smi=1
<a name="l00437"></a>00437             aet = smi*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a6e3b8c0a4f36c24d8af1238d81c33945">pet_cur</a>
<a name="l00438"></a>00438             Bvalue = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a05f361eeb6230c24a1d7718395c71b3b">rain_cur</a>-aet
<a name="l00439"></a>00439         
<a name="l00440"></a>00440             K = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a57edd92079beb2139337a15db7868836">theta2kr</a>(theta,thetar,thetas,m,l,Ks)
<a name="l00441"></a>00441             smc = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#aa83c63285e5f5efad117686fc9683af5">smcf</a>(theta,thetar,thetas,alpha,m,n)
<a name="l00442"></a>00442             psi = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a0445a3a4323788bfeb39b0dd4bf22e4a">theta2psi</a>(theta,thetar,thetas,m,n,alpha)
<a name="l00443"></a>00443                         
<a name="l00444"></a>00444             <span class="comment">#flux boundary condition at the top</span>
<a name="l00445"></a>00445             Kmid = np.empty(nz+1)        
<a name="l00446"></a>00446             Kmid[0] = 0
<a name="l00447"></a>00447             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l00448"></a>00448                 Kmid[i] = 0.5*(K[i]+K[i-1])
<a name="l00449"></a>00449             Kmid[nz] = K[nz-1]
<a name="l00450"></a>00450             
<a name="l00451"></a>00451             <span class="comment">#Setting the coefficient for the internal nodes</span>
<a name="l00452"></a>00452             A = np.empty(nz)
<a name="l00453"></a>00453             B = np.empty(nz)
<a name="l00454"></a>00454             C = np.empty(nz)
<a name="l00455"></a>00455             D = np.empty(nz)
<a name="l00456"></a>00456             dz = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a54a1bf0c73e3c943a5aadf73751e058b">dz</a>
<a name="l00457"></a>00457             dz2 = dz**2
<a name="l00458"></a>00458             
<a name="l00459"></a>00459             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(nz):
<a name="l00460"></a>00460                 A[i] = -(Kmid[i]/dz2)
<a name="l00461"></a>00461                 B[i] = smc[i]/dt+(Kmid[i+1]+Kmid[i])/dz2
<a name="l00462"></a>00462                 C[i] = A[i]
<a name="l00463"></a>00463                 D[i] = smc[i]*psi[i]/dt-(Kmid[i+1]-Kmid[i])/dz
<a name="l00464"></a>00464             <span class="comment"># setting the coefficient for the top bc (flux boundary)</span>
<a name="l00465"></a>00465             i = 0        
<a name="l00466"></a>00466             A[0] = 0
<a name="l00467"></a>00467             B[0] = smc[i]/dt+(Kmid[1])/dz2
<a name="l00468"></a>00468             D[0] = smc[i]*psi[i]/dt+(Bvalue-Kmid[1])/dz
<a name="l00469"></a>00469             
<a name="l00470"></a>00470             <span class="comment"># setting the coefficient for the bottom bc: gravity drainage</span>
<a name="l00471"></a>00471             B[nz-1] = smc[nz-1]/dt+(Kmid[nz])/dz2
<a name="l00472"></a>00472             C[nz-1] = 0
<a name="l00473"></a>00473             D[nz-1] = smc[nz-1]*psi[nz-1]/dt-(Kmid[nz]-Kmid[nz-1])/dz
<a name="l00474"></a>00474             
<a name="l00475"></a>00475             <span class="comment"># Solving using the thomas algorithm</span>
<a name="l00476"></a>00476             beta = np.empty(nz)
<a name="l00477"></a>00477             gamma = np.empty(nz)
<a name="l00478"></a>00478             u = np.empty(nz)
<a name="l00479"></a>00479             beta[0] = B[0]
<a name="l00480"></a>00480             gamma[0] = D[0]/beta[0]
<a name="l00481"></a>00481             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l00482"></a>00482                 beta[i] = B[i]-(A[i]*C[i-1])/(beta[i-1])
<a name="l00483"></a>00483                 gamma[i] = (D[i]-A[i]*gamma[i-1])/(beta[i])
<a name="l00484"></a>00484             
<a name="l00485"></a>00485             u[nz-1] = gamma[nz-1]
<a name="l00486"></a>00486             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(nz-2,-1,-1):
<a name="l00487"></a>00487                 u[i] = gamma[i]-(C[i]*u[i+1])/beta[i]
<a name="l00488"></a>00488             
<a name="l00489"></a>00489             <span class="comment"># flux computation between nodes</span>
<a name="l00490"></a>00490             J = np.empty(nz+1)
<a name="l00491"></a>00491             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l00492"></a>00492                 J[i] = Kmid[i]*(1-(u[i]-u[i-1])/dz)
<a name="l00493"></a>00493             J[0] = Bvalue
<a name="l00494"></a>00494             J[nz] = Kmid[nz]
<a name="l00495"></a>00495             
<a name="l00496"></a>00496             J[nz] = J[nz]
<a name="l00497"></a>00497             <span class="comment"># flux updating</span>
<a name="l00498"></a>00498             flux = np.diff(J)*dt/dz
<a name="l00499"></a>00499             theta = theta - flux
<a name="l00500"></a>00500             
<a name="l00501"></a>00501             <span class="keywordflow">if</span> theta[0]&gt;thetas:
<a name="l00502"></a>00502                 theta[theta&gt;thetas] = 0.99*thetas
<a name="l00503"></a>00503                         
<a name="l00504"></a>00504             aet_day += aet*dt 
<a name="l00505"></a>00505             recharge_day += J[nz]*dt
<a name="l00506"></a>00506                             
<a name="l00507"></a>00507         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ade9c77173908940774a8581b18cd7a74">theta</a> = theta        
<a name="l00508"></a>00508                       
<a name="l00509"></a>00509         <span class="comment"># write the output</span>
<a name="l00510"></a>00510         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a60c00cb4c4e4c0934eec2b046b970e72">nc_year</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>] = (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a4473319b3c1ac79c3a1f65bb59334506">cur_year</a>)
<a name="l00511"></a>00511         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a8040d34cfb0fec5d81c37d8be3e0c569">nc_doy</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>] = (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a5dd83c846f2ef3fa7d7a3656e85df861">cur_doy</a>)
<a name="l00512"></a>00512         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#af7a9b280883bbb1ced757befda5155a3">nc_sm</a>[:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>+1] = theta
<a name="l00513"></a>00513         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#af46617063a50d73918d2f3b82163612e">nc_recharge</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>] = recharge_day
<a name="l00514"></a>00514         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#acdbad2c6ebdb951bf36b1177eb2d9023">nc_aet</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a>] = aet_day
<a name="l00515"></a>00515         
<a name="l00516"></a>00516         <span class="comment"># print progress</span>
<a name="l00517"></a>00517         <span class="keywordflow">if</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a> == int(0.25*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2702f64024162cdafa67eddfe7dc016d">max_t</a>):
<a name="l00518"></a>00518             output_message = <span class="stringliteral">&#39;25 % completed&#39;</span>
<a name="l00519"></a>00519             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ab183d7511d72f3612a15ef2e20bdfc37">_colored_output</a>(output_message, 32)
<a name="l00520"></a>00520         
<a name="l00521"></a>00521         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a> == int(0.5*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2702f64024162cdafa67eddfe7dc016d">max_t</a>):
<a name="l00522"></a>00522             output_message = <span class="stringliteral">&#39;50 % completed&#39;</span>
<a name="l00523"></a>00523             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ab183d7511d72f3612a15ef2e20bdfc37">_colored_output</a>(output_message, 32)
<a name="l00524"></a>00524         
<a name="l00525"></a>00525         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a> == int(0.75*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2702f64024162cdafa67eddfe7dc016d">max_t</a>):
<a name="l00526"></a>00526             output_message = <span class="stringliteral">&#39;75 % completed&#39;</span>
<a name="l00527"></a>00527             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ab183d7511d72f3612a15ef2e20bdfc37">_colored_output</a>(output_message, 32)
<a name="l00528"></a>00528         
<a name="l00529"></a>00529         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a07e964f9716167a0666c4613aa5d3f0f" title="run the model ########################">t</a> == self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#a2702f64024162cdafa67eddfe7dc016d">max_t</a>-1:
<a name="l00530"></a>00530             output_message = <span class="stringliteral">&#39;100 % completed&#39;</span>
<a name="l00531"></a>00531             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html#ab183d7511d72f3612a15ef2e20bdfc37">_colored_output</a>(output_message, 32)
<a name="l00532"></a>00532         <span class="comment">#print self.t</span>
<a name="l00533"></a>00533 
<a name="l00534"></a>00534 
<a name="l00535"></a>00535 
<a name="l00536"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html">00536</a> <span class="keyword">class </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html">RICHARDS_1D_ENKF</a>(<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html">RICHARDS_1D</a>):
<a name="l00537"></a>00537     <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00538"></a>00538 <span class="stringliteral">    This is the main class of the Ensemble Kalman Filter (EnKF)</span>
<a name="l00539"></a>00539 <span class="stringliteral">    coupled with the one dimensional unsaturated model based on the </span>
<a name="l00540"></a>00540 <span class="stringliteral">    RICHARDS equation. The model is given in the class RICHARDS_1D.</span>
<a name="l00541"></a>00541 <span class="stringliteral">        </span>
<a name="l00542"></a>00542 <span class="stringliteral">    This will read the input data,</span>
<a name="l00543"></a>00543 <span class="stringliteral">    do the processing</span>
<a name="l00544"></a>00544 <span class="stringliteral">    and then write the output files</span>
<a name="l00545"></a>00545 <span class="stringliteral">    </span>
<a name="l00546"></a>00546 <span class="stringliteral">    &quot;&quot;&quot;</span>
<a name="l00547"></a>00547     
<a name="l00548"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a0ecdc4fe3795eb7c03a7a6a8da2ef2c1">00548</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a0ecdc4fe3795eb7c03a7a6a8da2ef2c1">__init__</a>(self,input_file):
<a name="l00549"></a>00549         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00550"></a>00550 <span class="stringliteral">        Input:</span>
<a name="l00551"></a>00551 <span class="stringliteral">            input_file: the file which contains all the information</span>
<a name="l00552"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">00552</a> <span class="stringliteral">            including forcing and parameters.</span>
<a name="l00553"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db">00553</a> <span class="stringliteral">        &quot;&quot;&quot;</span>      
<a name="l00554"></a>00554         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8f8a96b18ea802ee6a5a393e82e1ce52">input_file</a> = input_file
<a name="l00555"></a>00555         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a> = 10
<a name="l00556"></a>00556         <span class="comment"># read the input data</span>
<a name="l00557"></a>00557         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ae8ccc391434e64f597cd5c1bd3fd2bd8">_read_input</a>()
<a name="l00558"></a>00558         
<a name="l00559"></a>00559         <span class="comment"># initialize the variables and output file</span>
<a name="l00560"></a>00560         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a572c60532440a3e4910ee48ec33beb99">initialize</a>()
<a name="l00561"></a>00561         
<a name="l00562"></a>00562         <span class="comment">################ run the model ########################</span>
<a name="l00563"></a>00563         <span class="keywordflow">for</span> t <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a7bc8547c928fbe6248931112423126ac">max_t</a>):
<a name="l00564"></a>00564             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a> = t
<a name="l00565"></a>00565               
<a name="l00566"></a>00566             <span class="comment"># get forcing data at current time step        </span>
<a name="l00567"></a>00567             self._get_forcing()
<a name="l00568"></a>00568             
<a name="l00569"></a>00569             <span class="comment"># perturb the soil par ensemble</span>
<a name="l00570"></a>00570             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a6275a135bbe0c25912452380e3c8afdf">_perturb_soil_par_ens</a>()
<a name="l00571"></a>00571                         
<a name="l00572"></a>00572             <span class="comment"># call the unsat module with ensemble</span>
<a name="l00573"></a>00573             <span class="keywordflow">for</span> ens <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>):
<a name="l00574"></a>00574                 self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a22d47d63767a8512e39063ba8b6ebd29">ens</a> = ens
<a name="l00575"></a>00575                 
<a name="l00576"></a>00576                 self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ab1e24f77aa9c7fdc135d5974d75d2979">_unsat_ens</a>()
<a name="l00577"></a>00577                 
<a name="l00578"></a>00578             <span class="comment"># ensemble kalmfan filter</span>
<a name="l00579"></a>00579             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a1516736501e5ab117bfaa0387fc60fe7">_enkf_par_depth</a>()
<a name="l00580"></a>00580 
<a name="l00581"></a>00581             
<a name="l00582"></a>00582             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#acbc85012c408a546b820c2b06b987b91">_write_output</a>()                
<a name="l00583"></a>00583                 
<a name="l00584"></a>00584                 
<a name="l00585"></a>00585         self.nc_file.close() <span class="comment"># close the output file</span>
<a name="l00586"></a>00586 
<a name="l00587"></a>00587 
<a name="l00588"></a>00588     <span class="keyword">def </span>_enkf(self):
<a name="l00589"></a>00589         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00590"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">00590</a> <span class="stringliteral">        ensemble kalman filter</span>
<a name="l00591"></a>00591 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00592"></a>00592         <span class="comment"># compute the covariance matrix of the state+par</span>
<a name="l00593"></a>00593         x = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>
<a name="l00594"></a>00594         x_bar = np.tile(x.mean(axis=0),(10,1))
<a name="l00595"></a>00595         x_x_bar = x-x_bar
<a name="l00596"></a>00596         cov_xx = np.dot(x_x_bar.T,x_x_bar)        
<a name="l00597"></a>00597         
<a name="l00598"></a>00598         <span class="comment"># get the measurement of the ssm at the current time</span>
<a name="l00599"></a>00599         <span class="comment"># and generate its ensemble and compute its covariance matrix</span>
<a name="l00600"></a>00600         e = np.zeros((self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>, self.no_layer))
<a name="l00601"></a>00601         e[:,0] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>[:,0] - self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a658f5b3f6a3d84ca68f63d285f9aa930">meas_ssm</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>]
<a name="l00602"></a>00602         e = e + 0.02*np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,self.no_layer))
<a name="l00603"></a>00603         cov_ee = np.dot(e.T, e)  
<a name="l00604"></a>00604         
<a name="l00605"></a>00605         <span class="comment"># compute kalaman gain</span>
<a name="l00606"></a>00606         K = np.dot(cov_xx,np.linalg.pinv(cov_xx+cov_ee))
<a name="l00607"></a>00607         
<a name="l00608"></a>00608         <span class="comment"># update the measurment</span>
<a name="l00609"></a>00609         d = np.zeros(self.no_layer)
<a name="l00610"></a>00610         usm = np.zeros((self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,self.no_layer))
<a name="l00611"></a>00611              
<a name="l00612"></a>00612         <span class="keywordflow">for</span> ens <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>):
<a name="l00613"></a>00613             d[0] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a658f5b3f6a3d84ca68f63d285f9aa930">meas_ssm</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] - x[ens,0]
<a name="l00614"></a>00614             usm[ens,:] = (x[ens,:] + np.dot(K,d))
<a name="l00615"></a>00615                
<a name="l00616"></a>00616         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a> = usm
<a name="l00617"></a>00617     
<a name="l00618"></a>00618     <span class="keyword">def </span>_enkf_par(self):
<a name="l00619"></a>00619         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00620"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a16969c9e2b57f9bf832c091715c5c1fe">00620</a> <span class="stringliteral">        ensemble kalman filter</span>
<a name="l00621"></a>00621 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00622"></a>00622         <span class="comment"># make the state vector which contains the soil moisture at different </span>
<a name="l00623"></a>00623         <span class="comment">#depths and soil parameters</span>
<a name="l00624"></a>00624         x = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a> + 0.0*np.random.normal(size=self.no_layer)
<a name="l00625"></a>00625         thetar = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>]
<a name="l00626"></a>00626         thetas = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>]
<a name="l00627"></a>00627         alpha = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>]
<a name="l00628"></a>00628         n = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l00629"></a>00629         Ks = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>]
<a name="l00630"></a>00630         l = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;l&#39;</span>]
<a name="l00631"></a>00631         soil_par = (np.vstack([thetar, thetas, alpha, n, Ks, l])).T
<a name="l00632"></a>00632         X = np.hstack([x, soil_par])
<a name="l00633"></a>00633         
<a name="l00634"></a>00634         <span class="comment"># compute the covariance matrix of the state+par</span>
<a name="l00635"></a>00635         X_bar = np.tile(X.mean(axis=0),(10,1))
<a name="l00636"></a>00636         X_X_bar = X-X_bar
<a name="l00637"></a>00637         cov_XX = np.dot(X_X_bar.T,X_X_bar) + 1e-6*np.eye(self.no_layer+6)
<a name="l00638"></a>00638         cov_XX = 0.5*(cov_XX + cov_XX.T)
<a name="l00639"></a>00639         
<a name="l00640"></a>00640         <span class="comment"># get the measurement of the ssm at the current time</span>
<a name="l00641"></a>00641         <span class="comment"># and generate its ensemble and compute its covariance matrix</span>
<a name="l00642"></a>00642         e = np.zeros((self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>, self.no_layer+6))
<a name="l00643"></a>00643         e[:,0] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a658f5b3f6a3d84ca68f63d285f9aa930">meas_ssm</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>-1] - self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>[:,0].mean()
<a name="l00644"></a>00644         
<a name="l00645"></a>00645         v = 0.025*np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,self.no_layer+6))
<a name="l00646"></a>00646         v = v-np.tile(v.mean(axis=0),(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,1))
<a name="l00647"></a>00647         ev = e + v
<a name="l00648"></a>00648         cov_ee = np.dot(e.T, e) + 1e-6*np.eye(self.no_layer+6)
<a name="l00649"></a>00649         cov_ee = 0.5*(cov_ee + cov_ee.T)
<a name="l00650"></a>00650         
<a name="l00651"></a>00651         <span class="comment"># compute kalaman gain</span>
<a name="l00652"></a>00652         K = np.dot(cov_XX, np.linalg.pinv(cov_XX+cov_ee))
<a name="l00653"></a>00653         
<a name="l00654"></a>00654         <span class="comment"># update the measurment</span>
<a name="l00655"></a>00655         v = np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,))
<a name="l00656"></a>00656         v = v-v.mean()
<a name="l00657"></a>00657         e[:,0] = e[:,0]+0.0005*v
<a name="l00658"></a>00658         K = 0.5*(K + K.T)
<a name="l00659"></a>00659         usm_par = X + np.dot(K,e.T).T      
<a name="l00660"></a>00660         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a16969c9e2b57f9bf832c091715c5c1fe">usm_par</a> = usm_par
<a name="l00661"></a>00661         <span class="comment">#if self.t&lt;=5:</span>
<a name="l00662"></a>00662         <span class="comment">#    print self.t, usm_par[:,self.no_layer+4]/soil_par[:,4]</span>
<a name="l00663"></a>00663         <span class="comment"># check for the range of the updated ensemble</span>
<a name="l00664"></a>00664         <span class="comment"># soil moisture</span>
<a name="l00665"></a>00665         theta_ens = usm_par[:,:self.no_layer]
<a name="l00666"></a>00666         theta_ens[theta_ens&lt;0] = 0
<a name="l00667"></a>00667         theta_ens[theta_ens&gt;1] = 1
<a name="l00668"></a>00668         v = 0.01*np.random.normal(size=(theta_ens.shape))
<a name="l00669"></a>00669         v = v-np.tile(v.mean(axis=0),(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,1))
<a name="l00670"></a>00670         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a> = theta_ens + v
<a name="l00671"></a>00671         <span class="comment"># soil parameters</span>
<a name="l00672"></a>00672         thetar = usm_par[:,self.no_layer+0]
<a name="l00673"></a>00673         thetas = usm_par[:,self.no_layer+1]
<a name="l00674"></a>00674         alpha = usm_par[:,self.no_layer+2]
<a name="l00675"></a>00675         n = usm_par[:,self.no_layer+3]
<a name="l00676"></a>00676         Ks = usm_par[:,self.no_layer+4]
<a name="l00677"></a>00677         l = usm_par[:,self.no_layer+5]
<a name="l00678"></a>00678         
<a name="l00679"></a>00679         <span class="comment"># check the convergence limit</span>
<a name="l00680"></a>00680         <span class="comment">#thetar[thetar&gt;self.thetar_max] = 1.5*soil_par[:,0]</span>
<a name="l00681"></a>00681         <span class="comment">#thetas[thetas&gt;self.thetas_max] = 1.5*soil_par[:,1]</span>
<a name="l00682"></a>00682         <span class="comment">#alpha[alpha&gt;self.alpha_max] = 1.5*soil_par[:,2]</span>
<a name="l00683"></a>00683         <span class="comment">#n[n&gt;self.n_max] = 1.5*soil_par[:,3]</span>
<a name="l00684"></a>00684         <span class="comment">#Ks[Ks&gt;self.Ks_max] = 1.5*soil_par[:,4]</span>
<a name="l00685"></a>00685         <span class="comment">#l[l&gt;self.l_max] = 1.5*soil_par[:,5]</span>
<a name="l00686"></a>00686         
<a name="l00687"></a>00687         <span class="comment">#thetar[thetar&lt;self.thetar_min]  = 0.75*soil_par[:,0]</span>
<a name="l00688"></a>00688         <span class="comment">#thetas[thetas&lt;self.thetas_min]  = 0.75*soil_par[:,1]</span>
<a name="l00689"></a>00689         <span class="comment">#alpha[alpha&lt;self.alpha_min]     = 0.75*soil_par[:,2]</span>
<a name="l00690"></a>00690         <span class="comment">#n[n&lt;self.n_min]                 = 0.75*soil_par[:,3]</span>
<a name="l00691"></a>00691         <span class="comment">#Ks[Ks&lt;self.Ks_min]              = 0.75*soil_par[:,4]</span>
<a name="l00692"></a>00692         <span class="comment">#l[l&lt;self.l_min]                 = 0.75*soil_par[:,5]     </span>
<a name="l00693"></a>00693         
<a name="l00694"></a>00694         thetar[thetar&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>
<a name="l00695"></a>00695         thetas[thetas&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>
<a name="l00696"></a>00696         alpha[alpha&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>
<a name="l00697"></a>00697         n[n&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>
<a name="l00698"></a>00698         Ks[Ks&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>
<a name="l00699"></a>00699         l[l&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>
<a name="l00700"></a>00700         
<a name="l00701"></a>00701         thetar[thetar&lt;self.thetar_min]  = self.thetar_min
<a name="l00702"></a>00702         thetas[thetas&lt;self.thetas_min]  = self.thetas_min
<a name="l00703"></a>00703         alpha[alpha&lt;self.alpha_min]     = self.alpha_min
<a name="l00704"></a>00704         n[n&lt;self.n_min]                 = self.n_min
<a name="l00705"></a>00705         <span class="comment">#Ks[Ks&lt;self.Ks_min]              = self.Ks_min</span>
<a name="l00706"></a>00706         l[l&lt;self.l_min]                 = self.l_min
<a name="l00707"></a>00707         
<a name="l00708"></a>00708         
<a name="l00709"></a>00709         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>] = thetar
<a name="l00710"></a>00710         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>] = thetas
<a name="l00711"></a>00711         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>] = alpha
<a name="l00712"></a>00712         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;n&#39;</span>] = n
<a name="l00713"></a>00713         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>] = Ks
<a name="l00714"></a>00714         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;l&#39;</span>] = l
<a name="l00715"></a>00715         
<a name="l00716"></a>00716         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac17215a8eada0b319161dd1bb6e5bf5e">K</a> = K
<a name="l00717"></a>00717         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a3462463b46e53b86030b0ea4d011a02e">cov_ee</a> = cov_ee
<a name="l00718"></a>00718         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a1da0d41f4ee6481de3720465882dbc74">cov_XX</a> = cov_XX
<a name="l00719"></a>00719         
<a name="l00720"></a>00720     <span class="keyword">def </span>_enkf_par_depth(self):
<a name="l00721"></a>00721         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00722"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ae5b44f38de84318e6b4abb4762f21b7e">00722</a> <span class="stringliteral">        ensemble kalman filter</span>
<a name="l00723"></a>00723 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00724"></a>00724         <span class="comment"># make the state vector which contains the soil moisture at different </span>
<a name="l00725"></a>00725         <span class="comment">#depths and soil parameters</span>
<a name="l00726"></a>00726         x = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a> + 0.0*np.random.normal(size=self.no_layer)
<a name="l00727"></a>00727         thetar = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>]
<a name="l00728"></a>00728         thetas = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>]
<a name="l00729"></a>00729         alpha = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>]
<a name="l00730"></a>00730         n = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l00731"></a>00731         Ks = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>]
<a name="l00732"></a>00732         l = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;l&#39;</span>]
<a name="l00733"></a>00733         soil_par = (np.vstack([thetar, thetas, alpha, n, Ks, l])).T
<a name="l00734"></a>00734         X = np.hstack([x, soil_par])
<a name="l00735"></a>00735         
<a name="l00736"></a>00736         <span class="comment"># compute the covariance matrix of the state+par</span>
<a name="l00737"></a>00737         X_bar = np.tile(X.mean(axis=0),(10,1))
<a name="l00738"></a>00738         X_X_bar = X-X_bar
<a name="l00739"></a>00739         cov_XX = np.dot(X_X_bar.T,X_X_bar) + 1e-6*np.eye(self.no_layer+6)
<a name="l00740"></a>00740         cov_XX = 0.5*(cov_XX + cov_XX.T)
<a name="l00741"></a>00741         
<a name="l00742"></a>00742         <span class="comment"># get the measurement of the ssm at the current time</span>
<a name="l00743"></a>00743         <span class="comment"># and generate its ensemble and compute its covariance matrix</span>
<a name="l00744"></a>00744         e = np.zeros((self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>, self.no_layer+6))
<a name="l00745"></a>00745         ev = np.zeros((self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>, self.no_layer+6))
<a name="l00746"></a>00746         obs = np.zeros(self.no_layer)
<a name="l00747"></a>00747         obs =  self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ad9430e1593229aa2f45d4a470e95fe9a">a</a>+self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a746da5eecb8fa2fbcabbcd2a6e0e7cae">b</a>*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a658f5b3f6a3d84ca68f63d285f9aa930">meas_ssm</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>-1]
<a name="l00748"></a>00748         
<a name="l00749"></a>00749         <span class="comment">#for i in range(int(self.no_layer/2)):</span>
<a name="l00750"></a>00750         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(20):
<a name="l00751"></a>00751             e[:,i] = obs[i] - self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>[:,i].mean()
<a name="l00752"></a>00752         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ae5b44f38de84318e6b4abb4762f21b7e">e</a> = e            
<a name="l00753"></a>00753         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(self.no_layer+6):
<a name="l00754"></a>00754             <span class="keywordflow">if</span> i&lt;self.no_layer:
<a name="l00755"></a>00755                 ev[:,i] = e[:,i] + 0.005*(i**2+1)*np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>))
<a name="l00756"></a>00756             <span class="keywordflow">else</span>:
<a name="l00757"></a>00757                 ev[:,i] = e[:,i] + 0.000025*np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>))
<a name="l00758"></a>00758         
<a name="l00759"></a>00759         <span class="comment">#v = 0.05*np.random.normal(size=(self.n_ens,self.no_layer+6))</span>
<a name="l00760"></a>00760         <span class="comment">#v = v-np.tile(v.mean(axis=0),(self.n_ens,1))</span>
<a name="l00761"></a>00761         <span class="comment">#ev = e + v</span>
<a name="l00762"></a>00762         cov_ee = np.dot(ev.T, ev) + 1e-6*np.eye(self.no_layer+6)
<a name="l00763"></a>00763         cov_ee = 0.5*(cov_ee + cov_ee.T)
<a name="l00764"></a>00764         
<a name="l00765"></a>00765         <span class="comment"># compute kalaman gain</span>
<a name="l00766"></a>00766         K = np.dot(cov_XX, np.linalg.pinv(cov_XX+cov_ee))
<a name="l00767"></a>00767         
<a name="l00768"></a>00768         <span class="comment"># update the measurment</span>
<a name="l00769"></a>00769         v = np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,))
<a name="l00770"></a>00770         v = v-v.mean()
<a name="l00771"></a>00771         e[:,0] = e[:,0]+0.0005*v
<a name="l00772"></a>00772         K = 0.5*(K + K.T)
<a name="l00773"></a>00773         usm_par = X + np.dot(K,e.T).T      
<a name="l00774"></a>00774         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a16969c9e2b57f9bf832c091715c5c1fe">usm_par</a> = usm_par
<a name="l00775"></a>00775         <span class="comment">#if self.t&lt;=5:</span>
<a name="l00776"></a>00776         <span class="comment">#    print self.t, usm_par[:,self.no_layer+4]/soil_par[:,4]</span>
<a name="l00777"></a>00777         <span class="comment"># check for the range of the updated ensemble</span>
<a name="l00778"></a>00778         <span class="comment"># soil moisture</span>
<a name="l00779"></a>00779         theta_ens = usm_par[:,:self.no_layer]
<a name="l00780"></a>00780         theta_ens[theta_ens&lt;0] = 0
<a name="l00781"></a>00781         theta_ens[theta_ens&gt;1] = 1
<a name="l00782"></a>00782         v = 0.01*np.random.normal(size=(theta_ens.shape))
<a name="l00783"></a>00783         v = v-np.tile(v.mean(axis=0),(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,1))
<a name="l00784"></a>00784         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a> = theta_ens + v
<a name="l00785"></a>00785         <span class="comment"># soil parameters</span>
<a name="l00786"></a>00786         thetar = usm_par[:,self.no_layer+0]
<a name="l00787"></a>00787         thetas = usm_par[:,self.no_layer+1]
<a name="l00788"></a>00788         alpha = usm_par[:,self.no_layer+2]
<a name="l00789"></a>00789         n = usm_par[:,self.no_layer+3]
<a name="l00790"></a>00790         Ks = usm_par[:,self.no_layer+4]
<a name="l00791"></a>00791         l = usm_par[:,self.no_layer+5]
<a name="l00792"></a>00792         
<a name="l00793"></a>00793         <span class="comment"># check the convergence limit</span>
<a name="l00794"></a>00794         <span class="comment">#thetar[thetar&gt;self.thetar_max] = 1.5*soil_par[:,0]</span>
<a name="l00795"></a>00795         <span class="comment">#thetas[thetas&gt;self.thetas_max] = 1.5*soil_par[:,1]</span>
<a name="l00796"></a>00796         <span class="comment">#alpha[alpha&gt;self.alpha_max] = 1.5*soil_par[:,2]</span>
<a name="l00797"></a>00797         <span class="comment">#n[n&gt;self.n_max] = 1.5*soil_par[:,3]</span>
<a name="l00798"></a>00798         <span class="comment">#Ks[Ks&gt;self.Ks_max] = 1.5*soil_par[:,4]</span>
<a name="l00799"></a>00799         <span class="comment">#l[l&gt;self.l_max] = 1.5*soil_par[:,5]</span>
<a name="l00800"></a>00800         
<a name="l00801"></a>00801         <span class="comment">#thetar[thetar&lt;self.thetar_min]  = 0.75*soil_par[:,0]</span>
<a name="l00802"></a>00802         <span class="comment">#thetas[thetas&lt;self.thetas_min]  = 0.75*soil_par[:,1]</span>
<a name="l00803"></a>00803         <span class="comment">#alpha[alpha&lt;self.alpha_min]     = 0.75*soil_par[:,2]</span>
<a name="l00804"></a>00804         <span class="comment">#n[n&lt;self.n_min]                 = 0.75*soil_par[:,3]</span>
<a name="l00805"></a>00805         <span class="comment">#Ks[Ks&lt;self.Ks_min]              = 0.75*soil_par[:,4]</span>
<a name="l00806"></a>00806         <span class="comment">#l[l&lt;self.l_min]                 = 0.75*soil_par[:,5]     </span>
<a name="l00807"></a>00807         
<a name="l00808"></a>00808         thetar[thetar&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>
<a name="l00809"></a>00809         thetas[thetas&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>
<a name="l00810"></a>00810         alpha[alpha&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>
<a name="l00811"></a>00811         n[n&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>
<a name="l00812"></a>00812         Ks[Ks&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>
<a name="l00813"></a>00813         l[l&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>
<a name="l00814"></a>00814         
<a name="l00815"></a>00815         thetar[thetar&lt;self.thetar_min]  = self.thetar_min
<a name="l00816"></a>00816         thetas[thetas&lt;self.thetas_min]  = self.thetas_min
<a name="l00817"></a>00817         alpha[alpha&lt;self.alpha_min]     = self.alpha_min
<a name="l00818"></a>00818         n[n&lt;self.n_min]                 = self.n_min
<a name="l00819"></a>00819         Ks[Ks&lt;self.Ks_min]              = self.Ks_min
<a name="l00820"></a>00820         l[l&lt;self.l_min]                 = self.l_min
<a name="l00821"></a>00821         
<a name="l00822"></a>00822         
<a name="l00823"></a>00823         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>] = thetar
<a name="l00824"></a>00824         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>] = thetas
<a name="l00825"></a>00825         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>] = alpha
<a name="l00826"></a>00826         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;n&#39;</span>] = n
<a name="l00827"></a>00827         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>] = Ks
<a name="l00828"></a>00828         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;l&#39;</span>] = l
<a name="l00829"></a>00829         
<a name="l00830"></a>00830         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac17215a8eada0b319161dd1bb6e5bf5e">K</a> = K
<a name="l00831"></a>00831         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a3462463b46e53b86030b0ea4d011a02e">cov_ee</a> = cov_ee
<a name="l00832"></a>00832         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a1da0d41f4ee6481de3720465882dbc74">cov_XX</a> = cov_XX
<a name="l00833"></a>00833 
<a name="l00834"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a572c60532440a3e4910ee48ec33beb99">00834</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a572c60532440a3e4910ee48ec33beb99">initialize</a>(self):
<a name="l00835"></a>00835         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00836"></a>00836 <span class="stringliteral">        this initializes all the required variables</span>
<a name="l00837"></a>00837 <span class="stringliteral">        and open the netcdf file for writting</span>
<a name="l00838"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a677ec7938edf1d2dcd429d8df5fbf0b4">00838</a> <span class="stringliteral">        also generates the initial ensemble of the soil hydraulic parameters</span>
<a name="l00839"></a>00839 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00840"></a>00840         max_t = int(self.final_time/self.dt_flux)
<a name="l00841"></a>00841         <span class="comment">#max_t = 56</span>
<a name="l00842"></a>00842         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a7bc8547c928fbe6248931112423126ac">max_t</a> = max_t
<a name="l00843"></a>00843         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac9ed7496a5b32d3e4e461360818504c9">iter_dt</a> = 1
<a name="l00844"></a>00844                         
<a name="l00845"></a>00845         <span class="comment"># open file for writing</span>
<a name="l00846"></a>00846         file = nc.NetCDFFile(self.ofile_name, <span class="stringliteral">&#39;w&#39;</span>)
<a name="l00847"></a>00847         setattr(file, <span class="stringliteral">&#39;title&#39;</span>, <span class="stringliteral">&#39;output of the model ambhas.richards&#39;</span>)
<a name="l00848"></a>00848         now = datetime.datetime.now()
<a name="l00849"></a>00849         setattr(file, <span class="stringliteral">&#39;description&#39;</span>, <span class="stringliteral">&#39;The model was run at %s&#39;</span>%(now.ctime()))
<a name="l00850"></a>00850         file.createDimension(<span class="stringliteral">&#39;depth&#39;</span>, self.no_layer)
<a name="l00851"></a>00851         file.createDimension(<span class="stringliteral">&#39;time&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a7bc8547c928fbe6248931112423126ac">max_t</a>+1)
<a name="l00852"></a>00852         file.createDimension(<span class="stringliteral">&#39;ensemble&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>)
<a name="l00853"></a>00853         
<a name="l00854"></a>00854         <span class="comment"># depth</span>
<a name="l00855"></a>00855         varDims = <span class="stringliteral">&#39;depth&#39;</span>,
<a name="l00856"></a>00856         depth = file.createVariable(<span class="stringliteral">&#39;depth&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00857"></a>00857         depth.units = <span class="stringliteral">&#39;m&#39;</span>
<a name="l00858"></a>00858         depth[:] = np.tile(self.dz,self.no_layer).cumsum()-self.dz/2
<a name="l00859"></a>00859         
<a name="l00860"></a>00860         <span class="comment"># time (year and doy)</span>
<a name="l00861"></a>00861         varDims = <span class="stringliteral">&#39;time&#39;</span>,
<a name="l00862"></a>00862         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a677ec7938edf1d2dcd429d8df5fbf0b4">nc_year</a> = file.createVariable(<span class="stringliteral">&#39;year&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00863"></a>00863         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a6893099ef2198bd53c903bc6022f1cb7">nc_doy</a> = file.createVariable(<span class="stringliteral">&#39;doy&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00864"></a>00864         
<a name="l00865"></a>00865         <span class="comment"># soil moisture</span>
<a name="l00866"></a>00866         varDims = <span class="stringliteral">&#39;ensemble&#39;</span>, <span class="stringliteral">&#39;depth&#39;</span>, <span class="stringliteral">&#39;time&#39;</span>
<a name="l00867"></a>00867         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a97dd2f3ec9dfa80897c130c41cb32045">nc_sm</a> = file.createVariable(<span class="stringliteral">&#39;sm&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l00868"></a>00868         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a97dd2f3ec9dfa80897c130c41cb32045">nc_sm</a>.units = <span class="stringliteral">&#39;v/v&#39;</span>
<a name="l00869"></a>00869         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a97dd2f3ec9dfa80897c130c41cb32045">nc_sm</a>[:,:,0] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>
<a name="l00870"></a>00870         
<a name="l00871"></a>00871         <span class="comment"># recharge and aet</span>
<a name="l00872"></a>00872         varDims = <span class="stringliteral">&#39;time&#39;</span>,
<a name="l00873"></a>00873         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac01a9b6f5ed934b8cecf916debe65f57">nc_aet</a> = file.createVariable(<span class="stringliteral">&#39;aet&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00874"></a>00874         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac01a9b6f5ed934b8cecf916debe65f57">nc_aet</a>.units = <span class="stringliteral">&#39;mm&#39;</span>
<a name="l00875"></a>00875         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a27e53d39acde75790d361d073b9f19f2">nc_recharge</a> = file.createVariable(<span class="stringliteral">&#39;recharge&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00876"></a>00876         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a27e53d39acde75790d361d073b9f19f2">nc_recharge</a>.units = <span class="stringliteral">&#39;mm&#39;</span>
<a name="l00877"></a>00877 
<a name="l00878"></a>00878         <span class="comment"># recharge and aet</span>
<a name="l00879"></a>00879         varDims = <span class="stringliteral">&#39;ensemble&#39;</span>,<span class="stringliteral">&#39;time&#39;</span>
<a name="l00880"></a>00880         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ad711bb5c2231c41ea72d95c5b89eb0fc">nc_thetar</a> = file.createVariable(<span class="stringliteral">&#39;thetar&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00881"></a>00881         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ad711bb5c2231c41ea72d95c5b89eb0fc">nc_thetar</a>.units = <span class="stringliteral">&#39;v/v&#39;</span>
<a name="l00882"></a>00882         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a2b641d6fb1225d09dae3aaeec78b1eb9">nc_thetas</a> = file.createVariable(<span class="stringliteral">&#39;thetas&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00883"></a>00883         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a2b641d6fb1225d09dae3aaeec78b1eb9">nc_thetas</a>.units = <span class="stringliteral">&#39;v/v&#39;</span>    
<a name="l00884"></a>00884         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a6abc005b462e955071eb92e84e7e1abd">nc_alpha</a> = file.createVariable(<span class="stringliteral">&#39;alpha&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00885"></a>00885         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a6abc005b462e955071eb92e84e7e1abd">nc_alpha</a>.units = <span class="stringliteral">&#39;1/m&#39;</span>
<a name="l00886"></a>00886         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afbd4f9e7161f3ec5b490d5263ac7cc76">nc_n</a> = file.createVariable(<span class="stringliteral">&#39;n&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00887"></a>00887         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afbd4f9e7161f3ec5b490d5263ac7cc76">nc_n</a>.units = <span class="stringliteral">&#39;-&#39;</span> 
<a name="l00888"></a>00888         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4f1effc21e064b2e7f70cb818096ef8c">nc_Ks</a> = file.createVariable(<span class="stringliteral">&#39;Ks&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00889"></a>00889         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4f1effc21e064b2e7f70cb818096ef8c">nc_Ks</a>.units = <span class="stringliteral">&#39;m/s&#39;</span>
<a name="l00890"></a>00890         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa073385ed83280cf1d9b91e338187d7b">nc_l</a> = file.createVariable(<span class="stringliteral">&#39;l&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l00891"></a>00891         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa073385ed83280cf1d9b91e338187d7b">nc_l</a>.units = <span class="stringliteral">&#39;-&#39;</span> 
<a name="l00892"></a>00892         
<a name="l00893"></a>00893         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a5b06f90bfc9dade1717a1bfd07f184ff">nc_file</a> = file
<a name="l00894"></a>00894         
<a name="l00895"></a>00895         <span class="comment"># generate soil hydraulic parameters</span>
<a name="l00896"></a>00896         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aba360d7c58b58d648a2c6e8d06a61920">_generate_soil_par_ens</a>()
<a name="l00897"></a>00897 
<a name="l00898"></a>00898     <span class="keyword">def </span>_read_input(self):
<a name="l00899"></a>00899         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00900"></a>00900 <span class="stringliteral">        This checks if all the required input sheets are present in the xls file,</span>
<a name="l00901"></a>00901 <span class="stringliteral">        read the data from input file, which can be used later in other functions</span>
<a name="l00902"></a>00902 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00903"></a>00903         
<a name="l00904"></a>00904         <span class="comment"># list of required files in the input directory</span>
<a name="l00905"></a>00905         input_sheets = [<span class="stringliteral">&#39;ind&#39;</span>, <span class="stringliteral">&#39;forcing&#39;</span>, <span class="stringliteral">&#39;initial_condition&#39;</span>, <span class="stringliteral">&#39;units&#39;</span>, <span class="stringliteral">&#39;temporal_info&#39;</span>,
<a name="l00906"></a>00906                        <span class="stringliteral">&#39;spatial_info&#39;</span>, <span class="stringliteral">&#39;soil_hyd_par_ens&#39;</span>, <span class="stringliteral">&#39;output_par&#39;</span>]
<a name="l00907"></a>00907         
<a name="l00908"></a>00908         <span class="comment"># check if all the required sheets are present or not</span>
<a name="l00909"></a>00909         self._check_sheets(input_sheets, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8f8a96b18ea802ee6a5a393e82e1ce52">input_file</a>)
<a name="l00910"></a>00910         
<a name="l00911"></a>00911         <span class="comment"># read the legend</span>
<a name="l00912"></a>00912         self._read_ind()
<a name="l00913"></a>00913         
<a name="l00914"></a>00914         <span class="comment"># read the spatial data</span>
<a name="l00915"></a>00915         self._read_spatial()
<a name="l00916"></a>00916         
<a name="l00917"></a>00917         <span class="comment"># read the temporal data</span>
<a name="l00918"></a>00918         self._read_temporal()
<a name="l00919"></a>00919 
<a name="l00920"></a>00920         <span class="comment"># read the units </span>
<a name="l00921"></a>00921         self._read_units()
<a name="l00922"></a>00922         
<a name="l00923"></a>00923         <span class="comment"># read the initial condition</span>
<a name="l00924"></a>00924         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9f03fd23d9ebd50ea5dac530010aa014">_read_initial_condition</a>()
<a name="l00925"></a>00925         
<a name="l00926"></a>00926         <span class="comment"># read the ensemble information for the shp</span>
<a name="l00927"></a>00927         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a3969fdbf6997f994aa116056679bb66e">_read_shp_ens</a>()
<a name="l00928"></a>00928         
<a name="l00929"></a>00929         <span class="comment"># read the forcing infomation</span>
<a name="l00930"></a>00930         self._read_forcing()
<a name="l00931"></a>00931         
<a name="l00932"></a>00932         <span class="comment"># read the outfile name</span>
<a name="l00933"></a>00933         self._read_ofile_name()
<a name="l00934"></a>00934         
<a name="l00935"></a>00935         <span class="comment"># read the measured data</span>
<a name="l00936"></a>00936         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ab8c4ca3b305297d5414f03d4a8ddcd07">_read_measured</a>() 
<a name="l00937"></a>00937         
<a name="l00938"></a>00938         <span class="comment"># read a,b for the profile soil moisture relationship</span>
<a name="l00939"></a>00939         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa3a319618ad0677bdc895d6f93fbead4">_read_ab</a>()
<a name="l00940"></a>00940         
<a name="l00941"></a>00941         <span class="comment"># print the reading status</span>
<a name="l00942"></a>00942         output_message = <span class="stringliteral">&#39;Input data reading completed sucessfully&#39;</span>
<a name="l00943"></a>00943         self._colored_output(output_message, 32)
<a name="l00944"></a>00944     
<a name="l00945"></a>00945     <span class="keyword">def </span>_read_ab(self):
<a name="l00946"></a>00946         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00947"></a>00947 <span class="stringliteral">        read the intercept and slope of the relationship of the surface soil </span>
<a name="l00948"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a746da5eecb8fa2fbcabbcd2a6e0e7cae">00948</a> <span class="stringliteral">        moisture with the profile soil moisture</span>
<a name="l00949"></a>00949 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00950"></a>00950         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8f8a96b18ea802ee6a5a393e82e1ce52">input_file</a>)
<a name="l00951"></a>00951         sheet = book.sheet_by_name(<span class="stringliteral">&#39;ab&#39;</span>)
<a name="l00952"></a>00952         
<a name="l00953"></a>00953         data_len = sheet.nrows-1
<a name="l00954"></a>00954         a = np.zeros(data_len)
<a name="l00955"></a>00955         b = np.zeros(data_len)
<a name="l00956"></a>00956                
<a name="l00957"></a>00957         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> xrange(data_len):
<a name="l00958"></a>00958             a[i] = sheet.cell_value(i+1,1)
<a name="l00959"></a>00959             b[i] = sheet.cell_value(i+1,2)
<a name="l00960"></a>00960                     
<a name="l00961"></a>00961         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ad9430e1593229aa2f45d4a470e95fe9a">a</a> = a
<a name="l00962"></a>00962         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a746da5eecb8fa2fbcabbcd2a6e0e7cae">b</a> = b
<a name="l00963"></a>00963         
<a name="l00964"></a>00964     
<a name="l00965"></a>00965     <span class="keyword">def </span>_read_measured(self):
<a name="l00966"></a>00966         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00967"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a658f5b3f6a3d84ca68f63d285f9aa930">00967</a> <span class="stringliteral">        read the measured surface soil moisture (ssm) data</span>
<a name="l00968"></a>00968 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00969"></a>00969         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8f8a96b18ea802ee6a5a393e82e1ce52">input_file</a>)
<a name="l00970"></a>00970         sheet = book.sheet_by_name(<span class="stringliteral">&#39;forcing&#39;</span>)
<a name="l00971"></a>00971         
<a name="l00972"></a>00972         data_len = sheet.nrows-1
<a name="l00973"></a>00973         meas_ssm = np.zeros(data_len)
<a name="l00974"></a>00974         j = self.ind[<span class="stringliteral">&#39;meas_sm&#39;</span>]
<a name="l00975"></a>00975         
<a name="l00976"></a>00976         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> xrange(data_len):
<a name="l00977"></a>00977             meas_ssm[i] = sheet.cell_value(i+1,3+j)
<a name="l00978"></a>00978                     
<a name="l00979"></a>00979         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a658f5b3f6a3d84ca68f63d285f9aa930">meas_ssm</a> = meas_ssm
<a name="l00980"></a>00980     
<a name="l00981"></a>00981     <span class="keyword">def </span>_read_initial_condition(self):
<a name="l00982"></a>00982         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00983"></a>00983 <span class="stringliteral">        read initial condition</span>
<a name="l00984"></a>00984 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00985"></a>00985         <span class="comment">#get the row number from the ind</span>
<a name="l00986"></a>00986         j = self.ind[<span class="stringliteral">&#39;initial_condition&#39;</span>]
<a name="l00987"></a>00987         
<a name="l00988"></a>00988         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8f8a96b18ea802ee6a5a393e82e1ce52">input_file</a>)
<a name="l00989"></a>00989         sheet = book.sheet_by_name(<span class="stringliteral">&#39;initial_condition&#39;</span>)
<a name="l00990"></a>00990         theta_0 = sheet.cell_value(j,1)
<a name="l00991"></a>00991         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a> = theta_0 + 0.05*np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,self.no_layer))
<a name="l00992"></a>00992     
<a name="l00993"></a>00993     <span class="keyword">def </span>_read_shp_ens(self):
<a name="l00994"></a>00994         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00995"></a>00995 <span class="stringliteral">        read the information about the ensemble of the soil hydraulic parameters</span>
<a name="l00996"></a>00996 <span class="stringliteral">        the information being read is the min, max, best estimate and its uncertainty</span>
<a name="l00997"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">00997</a> <span class="stringliteral">        </span>
<a name="l00998"></a>00998 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l00999"></a>00999         <span class="comment">#get the row number from the ind</span>
<a name="l01000"></a>01000         j = self.ind[<span class="stringliteral">&#39;soil_hyd_par_ens&#39;</span>]
<a name="l01001"></a>01001         
<a name="l01002"></a>01002         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8f8a96b18ea802ee6a5a393e82e1ce52">input_file</a>)
<a name="l01003"></a>01003         sheet = book.sheet_by_name(<span class="stringliteral">&#39;soil_hyd_par_ens&#39;</span>)
<a name="l01004"></a>01004         self.thetar_min, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a> = sheet.cell_value(j+1,1), sheet.cell_value(j+1,7)
<a name="l01005"></a>01005         self.thetas_min, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a> = sheet.cell_value(j+1,2), sheet.cell_value(j+1,8)
<a name="l01006"></a>01006         self.alpha_min, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a> = sheet.cell_value(j+1,3), sheet.cell_value(j+1,9)
<a name="l01007"></a>01007         self.n_min, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a> = sheet.cell_value(j+1,4), sheet.cell_value(j+1,10)
<a name="l01008"></a>01008         self.Ks_min, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a> = sheet.cell_value(j+1,5), sheet.cell_value(j+1,11)
<a name="l01009"></a>01009         self.l_min, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a> = sheet.cell_value(j+1,6), sheet.cell_value(j+1,12)
<a name="l01010"></a>01010         
<a name="l01011"></a>01011         shp_ens = {}
<a name="l01012"></a>01012         shp_ens[<span class="stringliteral">&#39;thetar&#39;</span>] = sheet.cell_value(j+1,13), sheet.cell_value(j+1,19)
<a name="l01013"></a>01013         shp_ens[<span class="stringliteral">&#39;thetas&#39;</span>] = sheet.cell_value(j+1,14), sheet.cell_value(j+1,20)
<a name="l01014"></a>01014         shp_ens[<span class="stringliteral">&#39;alpha&#39;</span>]  = sheet.cell_value(j+1,15), sheet.cell_value(j+1,21)
<a name="l01015"></a>01015         shp_ens[<span class="stringliteral">&#39;n&#39;</span>]      = sheet.cell_value(j+1,16), sheet.cell_value(j+1,22)
<a name="l01016"></a>01016         shp_ens[<span class="stringliteral">&#39;Ks&#39;</span>]     = sheet.cell_value(j+1,17), sheet.cell_value(j+1,23)
<a name="l01017"></a>01017         shp_ens[<span class="stringliteral">&#39;l&#39;</span>]      = sheet.cell_value(j+1,18), sheet.cell_value(j+1,24)
<a name="l01018"></a>01018                 
<a name="l01019"></a>01019         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a> = shp_ens
<a name="l01020"></a>01020     
<a name="l01021"></a>01021     
<a name="l01022"></a>01022     <span class="keyword">def </span>_generate_soil_par_ens(self):
<a name="l01023"></a>01023         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01024"></a>01024 <span class="stringliteral">        this uses the LHS to generate the ensemble of the parameters</span>
<a name="l01025"></a>01025 <span class="stringliteral">        </span>
<a name="l01026"></a>01026 <span class="stringliteral">        this also computes the perturbation needed to perturb the parameters</span>
<a name="l01027"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">01027</a> <span class="stringliteral">        which is done in another function</span>
<a name="l01028"></a>01028 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01029"></a>01029                
<a name="l01030"></a>01030         <span class="comment">#gaussian perturbation</span>
<a name="l01031"></a>01031         v = np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,6))
<a name="l01032"></a>01032         v = v-np.tile(v.mean(axis=0),(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,1))
<a name="l01033"></a>01033         thetar = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>][0]  + v[:,0]*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>][1]
<a name="l01034"></a>01034         thetas = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>][0]  + v[:,1]*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>][1]
<a name="l01035"></a>01035         alpha = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>][0]    + v[:,2]*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>][1]
<a name="l01036"></a>01036         n = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;n&#39;</span>][0]            + v[:,3]*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;n&#39;</span>][1]
<a name="l01037"></a>01037         Ks = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>][0]          + v[:,4]*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>][1]
<a name="l01038"></a>01038         l = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;l&#39;</span>][0]            + v[:,5]*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a8500d0b699b3ad0ebacfdf74359a578f">shp_ens</a>[<span class="stringliteral">&#39;l&#39;</span>][1]
<a name="l01039"></a>01039 
<a name="l01040"></a>01040         <span class="comment"># check for the range of generated parameters</span>
<a name="l01041"></a>01041         thetar[thetar&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>
<a name="l01042"></a>01042         thetas[thetas&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>
<a name="l01043"></a>01043         alpha[alpha&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>
<a name="l01044"></a>01044         n[n&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>
<a name="l01045"></a>01045         Ks[Ks&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>
<a name="l01046"></a>01046         l[l&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>
<a name="l01047"></a>01047         
<a name="l01048"></a>01048         thetar[thetar&lt;self.thetar_min]  = self.thetar_min
<a name="l01049"></a>01049         thetas[thetas&lt;self.thetas_min]  = self.thetas_min
<a name="l01050"></a>01050         alpha[alpha&lt;self.alpha_min]     = self.alpha_min
<a name="l01051"></a>01051         n[n&lt;self.n_min]                 = self.n_min
<a name="l01052"></a>01052         Ks[Ks&lt;self.Ks_min]              = self.Ks_min
<a name="l01053"></a>01053         l[l&lt;self.l_min]                 = self.l_min
<a name="l01054"></a>01054         
<a name="l01055"></a>01055         soil_par_ens = {}
<a name="l01056"></a>01056         soil_par_ens[<span class="stringliteral">&#39;thetar&#39;</span>] = thetar
<a name="l01057"></a>01057         soil_par_ens[<span class="stringliteral">&#39;thetas&#39;</span>] = thetas
<a name="l01058"></a>01058         soil_par_ens[<span class="stringliteral">&#39;alpha&#39;</span>] = alpha
<a name="l01059"></a>01059         soil_par_ens[<span class="stringliteral">&#39;n&#39;</span>] = n
<a name="l01060"></a>01060         soil_par_ens[<span class="stringliteral">&#39;Ks&#39;</span>] = Ks
<a name="l01061"></a>01061         soil_par_ens[<span class="stringliteral">&#39;l&#39;</span>] = l
<a name="l01062"></a>01062         
<a name="l01063"></a>01063         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a> = soil_par_ens
<a name="l01064"></a>01064         
<a name="l01065"></a>01065         <span class="comment">#perturbation parameter</span>
<a name="l01066"></a>01066         soil_pert = {}
<a name="l01067"></a>01067         soil_pert[<span class="stringliteral">&#39;thetar&#39;</span>] = (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a> - self.thetar_min)*0.1/100.0
<a name="l01068"></a>01068         soil_pert[<span class="stringliteral">&#39;thetas&#39;</span>] = (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a> - self.thetas_min)*0.1/100.0
<a name="l01069"></a>01069         soil_pert[<span class="stringliteral">&#39;alpha&#39;</span>] =  (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a> - self.alpha_min)*0.1/100.0
<a name="l01070"></a>01070         soil_pert[<span class="stringliteral">&#39;n&#39;</span>] =      (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a> - self.n_min)*0.1/100.0
<a name="l01071"></a>01071         soil_pert[<span class="stringliteral">&#39;Ks&#39;</span>] =     (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a> - self.Ks_min)*0.1/100.0
<a name="l01072"></a>01072         soil_pert[<span class="stringliteral">&#39;l&#39;</span>] =      (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a> - self.l_min)*0.1/100.0
<a name="l01073"></a>01073         
<a name="l01074"></a>01074         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">soil_pert</a> = soil_pert
<a name="l01075"></a>01075 
<a name="l01076"></a>01076     <span class="keyword">def </span>_perturb_soil_par_ens(self):
<a name="l01077"></a>01077         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01078"></a>01078 <span class="stringliteral">        this functions perturb the soil hydraulic parameters </span>
<a name="l01079"></a>01079 <span class="stringliteral">        using the gaussian random variables</span>
<a name="l01080"></a>01080 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01081"></a>01081         v = np.random.normal(size=(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,6))
<a name="l01082"></a>01082         v = v-np.tile(v.mean(axis=0),(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4e38b7978e7b2e001f2fef63a2bb756d">n_ens</a>,1))
<a name="l01083"></a>01083         thetar = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>]+self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">soil_pert</a>[<span class="stringliteral">&#39;thetar&#39;</span>]*v[:,0]
<a name="l01084"></a>01084         thetas = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>]+self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">soil_pert</a>[<span class="stringliteral">&#39;thetas&#39;</span>]*v[:,1]
<a name="l01085"></a>01085         alpha = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>]+self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">soil_pert</a>[<span class="stringliteral">&#39;alpha&#39;</span>]*v[:,2]
<a name="l01086"></a>01086         n = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;n&#39;</span>]+self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">soil_pert</a>[<span class="stringliteral">&#39;n&#39;</span>]*v[:,3]
<a name="l01087"></a>01087         Ks = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>]+self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">soil_pert</a>[<span class="stringliteral">&#39;Ks&#39;</span>]*(v[:,4]-v[:,4].mean())
<a name="l01088"></a>01088         l = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;l&#39;</span>]+self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa9c56ab7b77e20ede15a75526a302ae1">soil_pert</a>[<span class="stringliteral">&#39;l&#39;</span>]*v[:,5]
<a name="l01089"></a>01089         
<a name="l01090"></a>01090 
<a name="l01091"></a>01091         <span class="comment"># check for the range of generated parameters</span>
<a name="l01092"></a>01092         thetar[thetar&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a925c18872e6998413f7cc8c6ad510a68">thetar_max</a>
<a name="l01093"></a>01093         thetas[thetas&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a060b417472eae60640654018b05a2d27">thetas_max</a>
<a name="l01094"></a>01094         alpha[alpha&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a28e172d19b032b47c29c3adf0ae57143">alpha_max</a>
<a name="l01095"></a>01095         n[n&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a81c29812703302fb10f4b346d026c4f9">n_max</a>
<a name="l01096"></a>01096         Ks[Ks&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a90ff64d1fc561985e394e3744e528a71">Ks_max</a>
<a name="l01097"></a>01097         l[l&gt;self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a290e7cb659097f6a1151ae7c3bcf8e99">l_max</a>
<a name="l01098"></a>01098         
<a name="l01099"></a>01099         thetar[thetar&lt;self.thetar_min]  = self.thetar_min
<a name="l01100"></a>01100         thetas[thetas&lt;self.thetas_min]  = self.thetas_min
<a name="l01101"></a>01101         alpha[alpha&lt;self.alpha_min]     = self.alpha_min
<a name="l01102"></a>01102         n[n&lt;self.n_min]                 = self.n_min
<a name="l01103"></a>01103         Ks[Ks&lt;self.Ks_min]              = self.Ks_min
<a name="l01104"></a>01104         l[l&lt;self.l_min]                 = self.l_min        
<a name="l01105"></a>01105 
<a name="l01106"></a>01106         soil_par_ens = {}
<a name="l01107"></a>01107         soil_par_ens[<span class="stringliteral">&#39;thetar&#39;</span>] = thetar
<a name="l01108"></a>01108         soil_par_ens[<span class="stringliteral">&#39;thetas&#39;</span>] = thetas
<a name="l01109"></a>01109         soil_par_ens[<span class="stringliteral">&#39;alpha&#39;</span>] = alpha
<a name="l01110"></a>01110         soil_par_ens[<span class="stringliteral">&#39;n&#39;</span>] = n
<a name="l01111"></a>01111         soil_par_ens[<span class="stringliteral">&#39;Ks&#39;</span>] = Ks
<a name="l01112"></a>01112         soil_par_ens[<span class="stringliteral">&#39;l&#39;</span>] = l        
<a name="l01113"></a>01113         
<a name="l01114"></a>01114         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a> = soil_par_ens
<a name="l01115"></a>01115 
<a name="l01116"></a>01116     
<a name="l01117"></a>01117     <span class="keyword">def </span>_unsat_ens(self):
<a name="l01118"></a>01118         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01119"></a>01119 <span class="stringliteral">        top boundary: atmoshpheric</span>
<a name="l01120"></a>01120 <span class="stringliteral">        bottom boundary: gravity drainage</span>
<a name="l01121"></a>01121 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01122"></a>01122         ens = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a22d47d63767a8512e39063ba8b6ebd29">ens</a>
<a name="l01123"></a>01123         thetar = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>][ens]
<a name="l01124"></a>01124         thetas = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>][ens]
<a name="l01125"></a>01125         alpha = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>][ens]
<a name="l01126"></a>01126         n = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;n&#39;</span>][ens]
<a name="l01127"></a>01127         l = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;l&#39;</span>][ens]
<a name="l01128"></a>01128         Ks = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>][ens]
<a name="l01129"></a>01129         m = 1-1/n
<a name="l01130"></a>01130         nz = self.no_layer
<a name="l01131"></a>01131         evap_0 = thetar+0.25*(thetas-thetar)
<a name="l01132"></a>01132         evap_1 = thetar+0.75*(thetas-thetar)
<a name="l01133"></a>01133         
<a name="l01134"></a>01134         <span class="comment"># read the soil moisture from the ensemble</span>
<a name="l01135"></a>01135         theta = 1.0*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>[ens]
<a name="l01136"></a>01136         
<a name="l01137"></a>01137         <span class="comment">#delta_theta = (np.abs(flux)).max()</span>
<a name="l01138"></a>01138         
<a name="l01139"></a>01139         iter_dt = max(24,int(np.ceil(self.rain_cur*self.dt_flux*1000/0.15)))
<a name="l01140"></a>01140         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac9ed7496a5b32d3e4e461360818504c9">iter_dt</a> = int(max(iter_dt,0.75*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac9ed7496a5b32d3e4e461360818504c9">iter_dt</a>))
<a name="l01141"></a>01141         
<a name="l01142"></a>01142         <span class="comment">#if self.t == 56:</span>
<a name="l01143"></a>01143         <span class="comment">#    self.iter_dt = int(self.iter_dt*6)</span>
<a name="l01144"></a>01144         <span class="comment">#print self.iter_dt</span>
<a name="l01145"></a>01145         
<a name="l01146"></a>01146         recharge_day = 0
<a name="l01147"></a>01147         aet_day = 0
<a name="l01148"></a>01148         
<a name="l01149"></a>01149         <span class="comment"># check for time step</span>
<a name="l01150"></a>01150         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac9ed7496a5b32d3e4e461360818504c9">iter_dt</a>):
<a name="l01151"></a>01151             dt = self.dt_flux/self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ac9ed7496a5b32d3e4e461360818504c9">iter_dt</a>
<a name="l01152"></a>01152             <span class="comment">#print dt</span>
<a name="l01153"></a>01153             <span class="comment"># top boundary value</span>
<a name="l01154"></a>01154             smi = (theta[0]-evap_0)/(evap_1-evap_0)
<a name="l01155"></a>01155             <span class="keywordflow">if</span> smi&lt;0: smi=0
<a name="l01156"></a>01156             <span class="keywordflow">if</span> smi&gt;1: smi=1
<a name="l01157"></a>01157             aet = smi*self.pet_cur
<a name="l01158"></a>01158             Bvalue = self.rain_cur-aet
<a name="l01159"></a>01159         
<a name="l01160"></a>01160             K = self.theta2kr(theta,thetar,thetas,m,l,Ks)
<a name="l01161"></a>01161             smc = self.smcf(theta,thetar,thetas,alpha,m,n)
<a name="l01162"></a>01162             psi = self.theta2psi(theta,thetar,thetas,m,n,alpha)
<a name="l01163"></a>01163             
<a name="l01164"></a>01164             
<a name="l01165"></a>01165             <span class="comment">#flux boundary condition at the top</span>
<a name="l01166"></a>01166             Kmid = np.empty(nz+1)        
<a name="l01167"></a>01167             Kmid[0] = 0
<a name="l01168"></a>01168             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l01169"></a>01169                 Kmid[i] = 0.5*(K[i]+K[i-1])
<a name="l01170"></a>01170             Kmid[nz] = K[nz-1]
<a name="l01171"></a>01171             
<a name="l01172"></a>01172             <span class="comment">#Setting the coefficient for the internal nodes</span>
<a name="l01173"></a>01173             A = np.empty(nz)
<a name="l01174"></a>01174             B = np.empty(nz)
<a name="l01175"></a>01175             C = np.empty(nz)
<a name="l01176"></a>01176             D = np.empty(nz)
<a name="l01177"></a>01177             dz = self.dz
<a name="l01178"></a>01178             dz2 = dz**2
<a name="l01179"></a>01179             
<a name="l01180"></a>01180             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(nz):
<a name="l01181"></a>01181                 A[i] = -(Kmid[i]/dz2)
<a name="l01182"></a>01182                 B[i] = smc[i]/dt+(Kmid[i+1]+Kmid[i])/dz2
<a name="l01183"></a>01183                 C[i] = A[i]
<a name="l01184"></a>01184                 D[i] = smc[i]*psi[i]/dt-(Kmid[i+1]-Kmid[i])/dz
<a name="l01185"></a>01185             <span class="comment"># setting the coefficient for the top bc (flux boundary)</span>
<a name="l01186"></a>01186             i = 0        
<a name="l01187"></a>01187             A[0] = 0
<a name="l01188"></a>01188             B[0] = smc[i]/dt+(Kmid[1])/dz2
<a name="l01189"></a>01189             D[0] = smc[i]*psi[i]/dt+(Bvalue-Kmid[1])/dz
<a name="l01190"></a>01190             
<a name="l01191"></a>01191             <span class="comment"># setting the coefficient for the bottom bc: gravity drainage</span>
<a name="l01192"></a>01192             B[nz-1] = smc[nz-1]/dt+(Kmid[nz])/dz2
<a name="l01193"></a>01193             C[nz-1] = 0
<a name="l01194"></a>01194             D[nz-1] = smc[nz-1]*psi[nz-1]/dt-(Kmid[nz]-Kmid[nz-1])/dz
<a name="l01195"></a>01195             
<a name="l01196"></a>01196             <span class="comment"># Solving using the thomas algorithm</span>
<a name="l01197"></a>01197             beta = np.empty(nz)
<a name="l01198"></a>01198             gamma = np.empty(nz)
<a name="l01199"></a>01199             u = np.empty(nz)
<a name="l01200"></a>01200             beta[0] = B[0]
<a name="l01201"></a>01201             gamma[0] = D[0]/beta[0]
<a name="l01202"></a>01202             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l01203"></a>01203                 beta[i] = B[i]-(A[i]*C[i-1])/(beta[i-1])
<a name="l01204"></a>01204                 gamma[i] = (D[i]-A[i]*gamma[i-1])/(beta[i])
<a name="l01205"></a>01205             
<a name="l01206"></a>01206             u[nz-1] = gamma[nz-1]
<a name="l01207"></a>01207             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(nz-2,-1,-1):
<a name="l01208"></a>01208                 u[i] = gamma[i]-(C[i]*u[i+1])/beta[i]
<a name="l01209"></a>01209             
<a name="l01210"></a>01210             <span class="comment"># flux computation between nodes</span>
<a name="l01211"></a>01211             J = np.empty(nz+1)
<a name="l01212"></a>01212             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l01213"></a>01213                 J[i] = Kmid[i]*(1-(u[i]-u[i-1])/dz)
<a name="l01214"></a>01214             J[0] = Bvalue
<a name="l01215"></a>01215             J[nz] = Kmid[nz]
<a name="l01216"></a>01216             
<a name="l01217"></a>01217             J[nz] = J[nz]
<a name="l01218"></a>01218             <span class="comment"># flux updating</span>
<a name="l01219"></a>01219             flux = np.diff(J)*dt/dz
<a name="l01220"></a>01220             theta = theta - flux
<a name="l01221"></a>01221             
<a name="l01222"></a>01222             theta[theta&gt;thetas] = 0.99*thetas
<a name="l01223"></a>01223             theta[theta&lt;thetar] = 1.01*thetar
<a name="l01224"></a>01224                         
<a name="l01225"></a>01225             aet_day += aet*dt 
<a name="l01226"></a>01226             recharge_day += J[nz]*dt           
<a name="l01227"></a>01227             
<a name="l01228"></a>01228                             
<a name="l01229"></a>01229         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>[ens] = theta        
<a name="l01230"></a>01230         
<a name="l01231"></a>01231     <span class="keyword">def </span>_write_output(self):
<a name="l01232"></a>01232         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01233"></a>01233 <span class="stringliteral">        this functions writes the output at each time step</span>
<a name="l01234"></a>01234 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01235"></a>01235         <span class="comment"># write the output</span>
<a name="l01236"></a>01236         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a677ec7938edf1d2dcd429d8df5fbf0b4">nc_year</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = (self.cur_year)
<a name="l01237"></a>01237         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a6893099ef2198bd53c903bc6022f1cb7">nc_doy</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = (self.cur_doy)
<a name="l01238"></a>01238         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a97dd2f3ec9dfa80897c130c41cb32045">nc_sm</a>[:,:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>+1] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afb0cd71a659acbdead58d7b478dc3122">theta_ens</a>
<a name="l01239"></a>01239         <span class="comment">#self.nc_recharge[self.t] = recharge_day</span>
<a name="l01240"></a>01240         <span class="comment">#self.nc_aet[self.t] = aet_day</span>
<a name="l01241"></a>01241         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#ad711bb5c2231c41ea72d95c5b89eb0fc">nc_thetar</a>[:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>]
<a name="l01242"></a>01242         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a2b641d6fb1225d09dae3aaeec78b1eb9">nc_thetas</a>[:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>]
<a name="l01243"></a>01243         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a6abc005b462e955071eb92e84e7e1abd">nc_alpha</a>[:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>]
<a name="l01244"></a>01244         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#afbd4f9e7161f3ec5b490d5263ac7cc76">nc_n</a>[:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l01245"></a>01245         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a4f1effc21e064b2e7f70cb818096ef8c">nc_Ks</a>[:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>]
<a name="l01246"></a>01246         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#aa073385ed83280cf1d9b91e338187d7b">nc_l</a>[:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a9ffb1a9fc95abf49ee53e0c88c2fe0db" title="run the model ########################">t</a>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__ENKF.html#a23474080a401ecf336487dbf7963bbfa">soil_par_ens</a>[<span class="stringliteral">&#39;l&#39;</span>]
<a name="l01247"></a>01247 
<a name="l01248"></a>01248 
<a name="l01249"></a>01249 
<a name="l01250"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html">01250</a> <span class="keyword">class </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html">RICHARDS_1D_GLUE</a>(<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D.html">RICHARDS_1D</a>):
<a name="l01251"></a>01251     <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01252"></a>01252 <span class="stringliteral">    This is the main class of the Ensemble Kalman Filter (EnKF)</span>
<a name="l01253"></a>01253 <span class="stringliteral">    coupled with the one dimensional unsaturated model based on the </span>
<a name="l01254"></a>01254 <span class="stringliteral">    RICHARDS equation. The model is given in the class RICHARDS_1D.</span>
<a name="l01255"></a>01255 <span class="stringliteral">        </span>
<a name="l01256"></a>01256 <span class="stringliteral">    This will read the input data,</span>
<a name="l01257"></a>01257 <span class="stringliteral">    do the processing</span>
<a name="l01258"></a>01258 <span class="stringliteral">    and then write the output files</span>
<a name="l01259"></a>01259 <span class="stringliteral">    </span>
<a name="l01260"></a>01260 <span class="stringliteral">    &quot;&quot;&quot;</span>
<a name="l01261"></a>01261     
<a name="l01262"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0db69bca7ee6bd0142105ffd18e30747">01262</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0db69bca7ee6bd0142105ffd18e30747">__init__</a>(self,input_file):
<a name="l01263"></a>01263         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01264"></a>01264 <span class="stringliteral">        Input:</span>
<a name="l01265"></a>01265 <span class="stringliteral">            input_file: the file which contains all the information</span>
<a name="l01266"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c59b9ae0ebac73e35397c84a004c8bd">01266</a> <span class="stringliteral">            including forcing and parameters.</span>
<a name="l01267"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305">01267</a> <span class="stringliteral">        &quot;&quot;&quot;</span>      
<a name="l01268"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0">01268</a>         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1c643e2120099a1854142e8dc647dc67">input_file</a> = input_file
<a name="l01269"></a>01269         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c59b9ae0ebac73e35397c84a004c8bd">n_ens</a> = 1000
<a name="l01270"></a>01270         
<a name="l01271"></a>01271         <span class="comment"># read the input data</span>
<a name="l01272"></a>01272         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a2bf14798e7e8b86843ce71a500534d06">_read_input</a>()
<a name="l01273"></a>01273         
<a name="l01274"></a>01274         <span class="comment"># initialize the variables and output file</span>
<a name="l01275"></a>01275         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0d7a55e5be5ca5c9c3e2ac9edc355784">initialize</a>()
<a name="l01276"></a>01276         
<a name="l01277"></a>01277         <span class="comment">########## run the GLUE ###########################</span>
<a name="l01278"></a>01278         <span class="keywordflow">for</span> ens <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c59b9ae0ebac73e35397c84a004c8bd">n_ens</a>):
<a name="l01279"></a>01279             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a> = ens
<a name="l01280"></a>01280             
<a name="l01281"></a>01281             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0a880cfc2a3255deb0ab22aec3f7542d">_shp_cur</a>()
<a name="l01282"></a>01282             
<a name="l01283"></a>01283             <span class="comment"># read the initial condition for ensemble everytime so that </span>
<a name="l01284"></a>01284             <span class="comment"># initial condition is same for all the ensemble</span>
<a name="l01285"></a>01285             self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1b60bd157c402e52d2653e121b746491">_read_initial_condition</a>()
<a name="l01286"></a>01286                         
<a name="l01287"></a>01287             <span class="comment">################ run the model ########################</span>
<a name="l01288"></a>01288             <span class="keywordflow">for</span> t <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a64367da04fa6c5ab496748cf7bbf621c">max_t</a>):
<a name="l01289"></a>01289                 self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a> = t
<a name="l01290"></a>01290                   
<a name="l01291"></a>01291                 <span class="comment"># get forcing data at current time step        </span>
<a name="l01292"></a>01292                 self._get_forcing()
<a name="l01293"></a>01293                 
<a name="l01294"></a>01294                 <span class="comment"># call the unsat module</span>
<a name="l01295"></a>01295                 self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a012c3e74f7b5d4cc5759ae08e05647bb">_unsat</a>()
<a name="l01296"></a>01296             
<a name="l01297"></a>01297             output_message = <span class="stringliteral">&#39;%d out of %d ensemble completed&#39;</span>%(ens,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c59b9ae0ebac73e35397c84a004c8bd">n_ens</a>)
<a name="l01298"></a>01298             self._colored_output(output_message, 41)
<a name="l01299"></a>01299             
<a name="l01300"></a>01300         self.nc_file.close() <span class="comment"># close the output file</span>
<a name="l01301"></a>01301 
<a name="l01302"></a>01302     <span class="keyword">def </span>_read_input(self):
<a name="l01303"></a>01303         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01304"></a>01304 <span class="stringliteral">        This checks if all the required input sheets are present in the xls file,</span>
<a name="l01305"></a>01305 <span class="stringliteral">        read the data from input file, which can be used later in other functions</span>
<a name="l01306"></a>01306 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01307"></a>01307         
<a name="l01308"></a>01308         <span class="comment"># list of required files in the input directory</span>
<a name="l01309"></a>01309         input_sheets = [<span class="stringliteral">&#39;ind&#39;</span>, <span class="stringliteral">&#39;forcing&#39;</span>, <span class="stringliteral">&#39;initial_condition&#39;</span>, <span class="stringliteral">&#39;units&#39;</span>, <span class="stringliteral">&#39;temporal_info&#39;</span>,
<a name="l01310"></a>01310                        <span class="stringliteral">&#39;spatial_info&#39;</span>, <span class="stringliteral">&#39;soil_hyd_par_ens&#39;</span>, <span class="stringliteral">&#39;output_par&#39;</span>]
<a name="l01311"></a>01311         
<a name="l01312"></a>01312         <span class="comment"># check if all the required sheets are present or not</span>
<a name="l01313"></a>01313         self._check_sheets(input_sheets, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1c643e2120099a1854142e8dc647dc67">input_file</a>)
<a name="l01314"></a>01314         
<a name="l01315"></a>01315         <span class="comment"># read the legend</span>
<a name="l01316"></a>01316         self._read_ind()
<a name="l01317"></a>01317         
<a name="l01318"></a>01318         <span class="comment"># read the spatial data</span>
<a name="l01319"></a>01319         self._read_spatial()
<a name="l01320"></a>01320         
<a name="l01321"></a>01321         <span class="comment"># read the temporal data</span>
<a name="l01322"></a>01322         self._read_temporal()
<a name="l01323"></a>01323 
<a name="l01324"></a>01324         <span class="comment"># read the units </span>
<a name="l01325"></a>01325         self._read_units()
<a name="l01326"></a>01326         
<a name="l01327"></a>01327         <span class="comment"># read the initial condition</span>
<a name="l01328"></a>01328         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1b60bd157c402e52d2653e121b746491">_read_initial_condition</a>()
<a name="l01329"></a>01329         
<a name="l01330"></a>01330         <span class="comment"># read the ensemble information for the shp</span>
<a name="l01331"></a>01331         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0300a305f3aa91f13133bda6a3f91755">_read_shp_ens</a>()
<a name="l01332"></a>01332         
<a name="l01333"></a>01333         <span class="comment"># read the forcing infomation</span>
<a name="l01334"></a>01334         self._read_forcing()
<a name="l01335"></a>01335         
<a name="l01336"></a>01336         <span class="comment"># read the outfile name</span>
<a name="l01337"></a>01337         self._read_ofile_name()
<a name="l01338"></a>01338         
<a name="l01339"></a>01339               
<a name="l01340"></a>01340                 
<a name="l01341"></a>01341         <span class="comment"># print the reading status</span>
<a name="l01342"></a>01342         output_message = <span class="stringliteral">&#39;Input data reading completed sucessfully&#39;</span>
<a name="l01343"></a>01343         self._colored_output(output_message, 32)
<a name="l01344"></a>01344     
<a name="l01345"></a>01345     
<a name="l01346"></a>01346     <span class="keyword">def </span>_read_shp_ens(self):
<a name="l01347"></a>01347         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01348"></a>01348 <span class="stringliteral">        read the information about the ensemble of the soil hydraulic parameters</span>
<a name="l01349"></a>01349 <span class="stringliteral">        the information being read is the min and max</span>
<a name="l01350"></a>01350 <span class="stringliteral">        </span>
<a name="l01351"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">01351</a> <span class="stringliteral">        this also uses LHS to generate the ensemble of the parameters        </span>
<a name="l01352"></a>01352 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01353"></a>01353         <span class="comment">#get the row number from the ind</span>
<a name="l01354"></a>01354         j = self.ind[<span class="stringliteral">&#39;soil_hyd_par_ens&#39;</span>]
<a name="l01355"></a>01355         
<a name="l01356"></a>01356         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1c643e2120099a1854142e8dc647dc67">input_file</a>)
<a name="l01357"></a>01357         sheet = book.sheet_by_name(<span class="stringliteral">&#39;soil_hyd_par_ens&#39;</span>)
<a name="l01358"></a>01358         thetar_min, thetar_max = sheet.cell_value(j+1,1), sheet.cell_value(j+1,7)
<a name="l01359"></a>01359         thetas_min, thetas_max = sheet.cell_value(j+1,2), sheet.cell_value(j+1,8)
<a name="l01360"></a>01360         alpha_min, alpha_max = sheet.cell_value(j+1,3), sheet.cell_value(j+1,9)
<a name="l01361"></a>01361         n_min, n_max = sheet.cell_value(j+1,4), sheet.cell_value(j+1,10)
<a name="l01362"></a>01362         Ks_min, Ks_max = sheet.cell_value(j+1,5), sheet.cell_value(j+1,11)
<a name="l01363"></a>01363         l_min, l_max = sheet.cell_value(j+1,6), sheet.cell_value(j+1,12)
<a name="l01364"></a>01364         
<a name="l01365"></a>01365         v = lhs(stats.uniform,[],(6,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c59b9ae0ebac73e35397c84a004c8bd">n_ens</a>))
<a name="l01366"></a>01366         
<a name="l01367"></a>01367         shp_ens = {}
<a name="l01368"></a>01368         shp_ens[<span class="stringliteral">&#39;thetar&#39;</span>] = thetar_min + (thetar_max-thetar_min)*v[0,:]
<a name="l01369"></a>01369         shp_ens[<span class="stringliteral">&#39;thetas&#39;</span>] = thetas_min + (thetas_max-thetas_min)*v[1,:]
<a name="l01370"></a>01370         shp_ens[<span class="stringliteral">&#39;alpha&#39;</span>]  = alpha_min + (alpha_max-alpha_min)*v[2,:]
<a name="l01371"></a>01371         shp_ens[<span class="stringliteral">&#39;n&#39;</span>]      = n_min + (n_max-n_min)*v[3,:]
<a name="l01372"></a>01372         shp_ens[<span class="stringliteral">&#39;Ks&#39;</span>]     = Ks_min + (Ks_max-Ks_min)*v[4,:]
<a name="l01373"></a>01373         shp_ens[<span class="stringliteral">&#39;l&#39;</span>]      = l_min + (l_max-l_min)*v[5,:]
<a name="l01374"></a>01374                 
<a name="l01375"></a>01375         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">shp_ens</a> = shp_ens
<a name="l01376"></a>01376     
<a name="l01377"></a>01377     <span class="keyword">def </span>_shp_cur(self):
<a name="l01378"></a>01378         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01379"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">01379</a> <span class="stringliteral">        read the current soil hydraulic parameters</span>
<a name="l01380"></a>01380 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01381"></a>01381         soil_par = {}
<a name="l01382"></a>01382         soil_par[<span class="stringliteral">&#39;thetar&#39;</span>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">shp_ens</a>[<span class="stringliteral">&#39;thetar&#39;</span>][self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>]           
<a name="l01383"></a>01383         soil_par[<span class="stringliteral">&#39;thetas&#39;</span>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">shp_ens</a>[<span class="stringliteral">&#39;thetas&#39;</span>][self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>]       
<a name="l01384"></a>01384         soil_par[<span class="stringliteral">&#39;alpha&#39;</span>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">shp_ens</a>[<span class="stringliteral">&#39;alpha&#39;</span>][self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>]       
<a name="l01385"></a>01385         soil_par[<span class="stringliteral">&#39;n&#39;</span>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">shp_ens</a>[<span class="stringliteral">&#39;n&#39;</span>][self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>]       
<a name="l01386"></a>01386         soil_par[<span class="stringliteral">&#39;Ks&#39;</span>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">shp_ens</a>[<span class="stringliteral">&#39;Ks&#39;</span>][self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>]       
<a name="l01387"></a>01387         soil_par[<span class="stringliteral">&#39;l&#39;</span>] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8e68be50b9c356c2246a39bf5e94b4af">shp_ens</a>[<span class="stringliteral">&#39;l&#39;</span>][self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>]  
<a name="l01388"></a>01388         soil_par[<span class="stringliteral">&#39;m&#39;</span>] = 1-1/soil_par[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l01389"></a>01389         
<a name="l01390"></a>01390         <span class="comment"># evaluate wilting point and field capacity</span>
<a name="l01391"></a>01391         soil_par[<span class="stringliteral">&#39;evap_1&#39;</span>] = self.psi2theta(-0.33, soil_par[<span class="stringliteral">&#39;thetar&#39;</span>], soil_par[<span class="stringliteral">&#39;thetas&#39;</span>], 
<a name="l01392"></a>01392                                soil_par[<span class="stringliteral">&#39;alpha&#39;</span>], soil_par[<span class="stringliteral">&#39;m&#39;</span>], soil_par[<span class="stringliteral">&#39;n&#39;</span>])
<a name="l01393"></a>01393         
<a name="l01394"></a>01394         soil_par[<span class="stringliteral">&#39;evap_0&#39;</span>] = self.psi2theta(-15, soil_par[<span class="stringliteral">&#39;thetar&#39;</span>], soil_par[<span class="stringliteral">&#39;thetas&#39;</span>], 
<a name="l01395"></a>01395                                soil_par[<span class="stringliteral">&#39;alpha&#39;</span>], soil_par[<span class="stringliteral">&#39;m&#39;</span>], soil_par[<span class="stringliteral">&#39;n&#39;</span>])
<a name="l01396"></a>01396         
<a name="l01397"></a>01397         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a> = soil_par
<a name="l01398"></a>01398         
<a name="l01399"></a>01399         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a46419513530d70c1c41dee95cff641dc">nc_thetar</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>] = soil_par[<span class="stringliteral">&#39;thetar&#39;</span>]
<a name="l01400"></a>01400         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a2cf6ded17b2b78ee83e3baf263a551ad">nc_thetas</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>] = soil_par[<span class="stringliteral">&#39;thetas&#39;</span>]
<a name="l01401"></a>01401         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a68b2d782921fe2dccac528081c9b71a0">nc_alpha</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>] = soil_par[<span class="stringliteral">&#39;alpha&#39;</span>]
<a name="l01402"></a>01402         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a50e39187a46d6cc8491761930f624f43">nc_n</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>] = soil_par[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l01403"></a>01403         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a7c736cecf880cc6fab74ea3990824906">nc_Ks</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>] = soil_par[<span class="stringliteral">&#39;Ks&#39;</span>]
<a name="l01404"></a>01404         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0473d96d27c986eab19a88ceaa3f734a">nc_l</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>] = soil_par[<span class="stringliteral">&#39;l&#39;</span>]
<a name="l01405"></a>01405                                
<a name="l01406"></a>01406     
<a name="l01407"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0d7a55e5be5ca5c9c3e2ac9edc355784">01407</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0d7a55e5be5ca5c9c3e2ac9edc355784">initialize</a>(self):
<a name="l01408"></a>01408         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01409"></a>01409 <span class="stringliteral">        this initializes all the required variables</span>
<a name="l01410"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a3abdd027570f1e8ee8b5a1c3937ebecd">01410</a> <span class="stringliteral">        and open the netcdf file for writting</span>
<a name="l01411"></a>01411 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01412"></a>01412         max_t = int(self.final_time/self.dt_flux)
<a name="l01413"></a>01413         
<a name="l01414"></a>01414         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a64367da04fa6c5ab496748cf7bbf621c">max_t</a> = max_t
<a name="l01415"></a>01415         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#ad3e8694827645301e38c703994194529">iter_dt</a> = 1
<a name="l01416"></a>01416                         
<a name="l01417"></a>01417         <span class="comment"># open file for writing</span>
<a name="l01418"></a>01418         file = nc.NetCDFFile(self.ofile_name, <span class="stringliteral">&#39;w&#39;</span>)
<a name="l01419"></a>01419         setattr(file, <span class="stringliteral">&#39;title&#39;</span>, <span class="stringliteral">&#39;output of the model ambhas.richards_glue&#39;</span>)
<a name="l01420"></a>01420         now = datetime.datetime.now()
<a name="l01421"></a>01421         setattr(file, <span class="stringliteral">&#39;description&#39;</span>, <span class="stringliteral">&#39;The model was run at %s&#39;</span>%(now.ctime()))
<a name="l01422"></a>01422         file.createDimension(<span class="stringliteral">&#39;depth&#39;</span>, self.no_layer)
<a name="l01423"></a>01423         file.createDimension(<span class="stringliteral">&#39;time&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a64367da04fa6c5ab496748cf7bbf621c">max_t</a>+1)
<a name="l01424"></a>01424         file.createDimension(<span class="stringliteral">&#39;ensemble&#39;</span>, self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c59b9ae0ebac73e35397c84a004c8bd">n_ens</a>)
<a name="l01425"></a>01425         
<a name="l01426"></a>01426         <span class="comment"># depth</span>
<a name="l01427"></a>01427         varDims = <span class="stringliteral">&#39;depth&#39;</span>,
<a name="l01428"></a>01428         depth = file.createVariable(<span class="stringliteral">&#39;depth&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l01429"></a>01429         depth.units = <span class="stringliteral">&#39;m&#39;</span>
<a name="l01430"></a>01430         depth[:] = np.tile(self.dz,self.no_layer).cumsum()-self.dz/2
<a name="l01431"></a>01431         
<a name="l01432"></a>01432         <span class="comment"># time (year and doy)</span>
<a name="l01433"></a>01433         varDims = <span class="stringliteral">&#39;time&#39;</span>,
<a name="l01434"></a>01434         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a3abdd027570f1e8ee8b5a1c3937ebecd">nc_year</a> = file.createVariable(<span class="stringliteral">&#39;year&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l01435"></a>01435         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#af00e43fdcd78dc1cf989e23c83576cf6">nc_doy</a> = file.createVariable(<span class="stringliteral">&#39;doy&#39;</span>, <span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l01436"></a>01436         
<a name="l01437"></a>01437         <span class="comment"># soil moisture</span>
<a name="l01438"></a>01438         varDims = <span class="stringliteral">&#39;ensemble&#39;</span>, <span class="stringliteral">&#39;depth&#39;</span>, <span class="stringliteral">&#39;time&#39;</span>
<a name="l01439"></a>01439         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a326ef99eeb5d1b35707081c5d556ccb7">nc_sm</a> = file.createVariable(<span class="stringliteral">&#39;sm&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>, varDims)
<a name="l01440"></a>01440         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a326ef99eeb5d1b35707081c5d556ccb7">nc_sm</a>.units = <span class="stringliteral">&#39;v/v&#39;</span>
<a name="l01441"></a>01441         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a326ef99eeb5d1b35707081c5d556ccb7">nc_sm</a>[:,:,0] = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1738f98e410c66a6e0ab4ec29bc0a264">theta</a>
<a name="l01442"></a>01442         
<a name="l01443"></a>01443         <span class="comment"># recharge and aet</span>
<a name="l01444"></a>01444         varDims = <span class="stringliteral">&#39;ensemble&#39;</span>,<span class="stringliteral">&#39;time&#39;</span>
<a name="l01445"></a>01445         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a919d83b479e4087a89d7216f8570b55b">nc_aet</a> = file.createVariable(<span class="stringliteral">&#39;aet&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01446"></a>01446         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a919d83b479e4087a89d7216f8570b55b">nc_aet</a>.units = <span class="stringliteral">&#39;mm&#39;</span>
<a name="l01447"></a>01447         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a45ec029f327b4811451d3db36b93c9a6">nc_recharge</a> = file.createVariable(<span class="stringliteral">&#39;recharge&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01448"></a>01448         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a45ec029f327b4811451d3db36b93c9a6">nc_recharge</a>.units = <span class="stringliteral">&#39;mm&#39;</span>
<a name="l01449"></a>01449         
<a name="l01450"></a>01450         <span class="comment"># soil_par</span>
<a name="l01451"></a>01451         varDims = <span class="stringliteral">&#39;ensemble&#39;</span>,
<a name="l01452"></a>01452         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a46419513530d70c1c41dee95cff641dc">nc_thetar</a> = file.createVariable(<span class="stringliteral">&#39;thetar&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01453"></a>01453         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a2cf6ded17b2b78ee83e3baf263a551ad">nc_thetas</a> = file.createVariable(<span class="stringliteral">&#39;thetas&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01454"></a>01454         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a68b2d782921fe2dccac528081c9b71a0">nc_alpha</a> = file.createVariable(<span class="stringliteral">&#39;alpha&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01455"></a>01455         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a50e39187a46d6cc8491761930f624f43">nc_n</a> = file.createVariable(<span class="stringliteral">&#39;n&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01456"></a>01456         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a7c736cecf880cc6fab74ea3990824906">nc_Ks</a> = file.createVariable(<span class="stringliteral">&#39;Ks&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01457"></a>01457         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a0473d96d27c986eab19a88ceaa3f734a">nc_l</a> = file.createVariable(<span class="stringliteral">&#39;l&#39;</span>,<span class="stringliteral">&#39;d&#39;</span>,varDims)
<a name="l01458"></a>01458         
<a name="l01459"></a>01459         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a114ef3c752202bbb85036518b0872a78">nc_file</a> = file
<a name="l01460"></a>01460     
<a name="l01461"></a>01461     <span class="keyword">def </span>_read_initial_condition(self):
<a name="l01462"></a>01462         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01463"></a><a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1738f98e410c66a6e0ab4ec29bc0a264">01463</a> <span class="stringliteral">        read initial condition</span>
<a name="l01464"></a>01464 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01465"></a>01465         <span class="comment">#get the row number from the ind</span>
<a name="l01466"></a>01466         j = self.ind[<span class="stringliteral">&#39;initial_condition&#39;</span>]
<a name="l01467"></a>01467         
<a name="l01468"></a>01468         book = xlrd.open_workbook(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1c643e2120099a1854142e8dc647dc67">input_file</a>)
<a name="l01469"></a>01469         sheet = book.sheet_by_name(<span class="stringliteral">&#39;initial_condition&#39;</span>)
<a name="l01470"></a>01470         
<a name="l01471"></a>01471         data_len = sheet.nrows-1
<a name="l01472"></a>01472         theta = np.zeros(data_len)
<a name="l01473"></a>01473         
<a name="l01474"></a>01474         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> xrange(data_len):
<a name="l01475"></a>01475             theta[i] = sheet.cell_value(i+1,j-1)
<a name="l01476"></a>01476         
<a name="l01477"></a>01477         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1738f98e410c66a6e0ab4ec29bc0a264">theta</a> = theta
<a name="l01478"></a>01478                 
<a name="l01479"></a>01479     
<a name="l01480"></a>01480     <span class="keyword">def </span>_unsat(self):
<a name="l01481"></a>01481         <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l01482"></a>01482 <span class="stringliteral">        top boundary: atmoshpheric</span>
<a name="l01483"></a>01483 <span class="stringliteral">        bottom boundary: gravity drainage</span>
<a name="l01484"></a>01484 <span class="stringliteral">        &quot;&quot;&quot;</span>
<a name="l01485"></a>01485                
<a name="l01486"></a>01486         thetar = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;thetar&#39;</span>]
<a name="l01487"></a>01487         thetas = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;thetas&#39;</span>]
<a name="l01488"></a>01488         alpha = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;alpha&#39;</span>]
<a name="l01489"></a>01489         n = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;n&#39;</span>]
<a name="l01490"></a>01490         m = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;m&#39;</span>]
<a name="l01491"></a>01491         l = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;l&#39;</span>]
<a name="l01492"></a>01492         Ks = self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;Ks&#39;</span>]
<a name="l01493"></a>01493         nz = self.no_layer
<a name="l01494"></a>01494                 
<a name="l01495"></a>01495         theta = 1.0*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1738f98e410c66a6e0ab4ec29bc0a264">theta</a>
<a name="l01496"></a>01496         
<a name="l01497"></a>01497         <span class="comment">#delta_theta = (np.abs(flux)).max()</span>
<a name="l01498"></a>01498         
<a name="l01499"></a>01499         iter_dt = max(24,int(np.ceil(self.rain_cur*self.dt_flux*1000/0.15)))
<a name="l01500"></a>01500         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#ad3e8694827645301e38c703994194529">iter_dt</a> = int(max(iter_dt,0.75*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#ad3e8694827645301e38c703994194529">iter_dt</a>))
<a name="l01501"></a>01501         
<a name="l01502"></a>01502         <span class="comment">#if self.t == 56:</span>
<a name="l01503"></a>01503         <span class="comment">#    self.iter_dt = int(self.iter_dt*6)</span>
<a name="l01504"></a>01504         <span class="comment">#print self.iter_dt</span>
<a name="l01505"></a>01505         
<a name="l01506"></a>01506         recharge_day = 0
<a name="l01507"></a>01507         aet_day = 0
<a name="l01508"></a>01508         
<a name="l01509"></a>01509         <span class="comment"># check for time step</span>
<a name="l01510"></a>01510         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#ad3e8694827645301e38c703994194529">iter_dt</a>):
<a name="l01511"></a>01511             dt = self.dt_flux/self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#ad3e8694827645301e38c703994194529">iter_dt</a>
<a name="l01512"></a>01512             <span class="comment">#print dt</span>
<a name="l01513"></a>01513             <span class="comment"># top boundary value</span>
<a name="l01514"></a>01514             smi = (self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1738f98e410c66a6e0ab4ec29bc0a264">theta</a>[0]-self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;evap_0&#39;</span>])/(self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;evap_1&#39;</span>]-self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a066e633a9d95cd646504440430124f57">soil_par</a>[<span class="stringliteral">&#39;evap_0&#39;</span>])
<a name="l01515"></a>01515             <span class="keywordflow">if</span> smi&lt;0: smi=0
<a name="l01516"></a>01516             <span class="keywordflow">if</span> smi&gt;1: smi=1
<a name="l01517"></a>01517             aet = smi*self.pet_cur
<a name="l01518"></a>01518             Bvalue = self.rain_cur-aet
<a name="l01519"></a>01519         
<a name="l01520"></a>01520             K = self.theta2kr(theta,thetar,thetas,m,l,Ks)
<a name="l01521"></a>01521             smc = self.smcf(theta,thetar,thetas,alpha,m,n)
<a name="l01522"></a>01522             psi = self.theta2psi(theta,thetar,thetas,m,n,alpha)
<a name="l01523"></a>01523                         
<a name="l01524"></a>01524             <span class="comment">#flux boundary condition at the top</span>
<a name="l01525"></a>01525             Kmid = np.empty(nz+1)        
<a name="l01526"></a>01526             Kmid[0] = 0
<a name="l01527"></a>01527             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l01528"></a>01528                 Kmid[i] = 0.5*(K[i]+K[i-1])
<a name="l01529"></a>01529             Kmid[nz] = K[nz-1]
<a name="l01530"></a>01530             
<a name="l01531"></a>01531             <span class="comment">#Setting the coefficient for the internal nodes</span>
<a name="l01532"></a>01532             A = np.empty(nz)
<a name="l01533"></a>01533             B = np.empty(nz)
<a name="l01534"></a>01534             C = np.empty(nz)
<a name="l01535"></a>01535             D = np.empty(nz)
<a name="l01536"></a>01536             dz = self.dz
<a name="l01537"></a>01537             dz2 = dz**2
<a name="l01538"></a>01538             
<a name="l01539"></a>01539             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(nz):
<a name="l01540"></a>01540                 A[i] = -(Kmid[i]/dz2)
<a name="l01541"></a>01541                 B[i] = smc[i]/dt+(Kmid[i+1]+Kmid[i])/dz2
<a name="l01542"></a>01542                 C[i] = A[i]
<a name="l01543"></a>01543                 D[i] = smc[i]*psi[i]/dt-(Kmid[i+1]-Kmid[i])/dz
<a name="l01544"></a>01544             <span class="comment"># setting the coefficient for the top bc (flux boundary)</span>
<a name="l01545"></a>01545             i = 0        
<a name="l01546"></a>01546             A[0] = 0
<a name="l01547"></a>01547             B[0] = smc[i]/dt+(Kmid[1])/dz2
<a name="l01548"></a>01548             D[0] = smc[i]*psi[i]/dt+(Bvalue-Kmid[1])/dz
<a name="l01549"></a>01549             
<a name="l01550"></a>01550             <span class="comment"># setting the coefficient for the bottom bc: gravity drainage</span>
<a name="l01551"></a>01551             B[nz-1] = smc[nz-1]/dt+(Kmid[nz])/dz2
<a name="l01552"></a>01552             C[nz-1] = 0
<a name="l01553"></a>01553             D[nz-1] = smc[nz-1]*psi[nz-1]/dt-(Kmid[nz]-Kmid[nz-1])/dz
<a name="l01554"></a>01554             
<a name="l01555"></a>01555             <span class="comment"># Solving using the thomas algorithm</span>
<a name="l01556"></a>01556             beta = np.empty(nz)
<a name="l01557"></a>01557             gamma = np.empty(nz)
<a name="l01558"></a>01558             u = np.empty(nz)
<a name="l01559"></a>01559             beta[0] = B[0]
<a name="l01560"></a>01560             gamma[0] = D[0]/beta[0]
<a name="l01561"></a>01561             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l01562"></a>01562                 beta[i] = B[i]-(A[i]*C[i-1])/(beta[i-1])
<a name="l01563"></a>01563                 gamma[i] = (D[i]-A[i]*gamma[i-1])/(beta[i])
<a name="l01564"></a>01564             
<a name="l01565"></a>01565             u[nz-1] = gamma[nz-1]
<a name="l01566"></a>01566             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(nz-2,-1,-1):
<a name="l01567"></a>01567                 u[i] = gamma[i]-(C[i]*u[i+1])/beta[i]
<a name="l01568"></a>01568             
<a name="l01569"></a>01569             <span class="comment"># flux computation between nodes</span>
<a name="l01570"></a>01570             J = np.empty(nz+1)
<a name="l01571"></a>01571             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(1,nz):
<a name="l01572"></a>01572                 J[i] = Kmid[i]*(1-(u[i]-u[i-1])/dz)
<a name="l01573"></a>01573             J[0] = Bvalue
<a name="l01574"></a>01574             J[nz] = Kmid[nz]
<a name="l01575"></a>01575             
<a name="l01576"></a>01576             J[nz] = J[nz]
<a name="l01577"></a>01577             <span class="comment"># flux updating</span>
<a name="l01578"></a>01578             flux = np.diff(J)*dt/dz
<a name="l01579"></a>01579             theta = theta - flux
<a name="l01580"></a>01580             
<a name="l01581"></a>01581             theta[theta&gt;thetas] = 0.99*thetas
<a name="l01582"></a>01582             theta[theta&lt;thetar] = 1.01*thetar
<a name="l01583"></a>01583                         
<a name="l01584"></a>01584             aet_day += aet*dt 
<a name="l01585"></a>01585             recharge_day += J[nz]*dt
<a name="l01586"></a>01586                             
<a name="l01587"></a>01587         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a1738f98e410c66a6e0ab4ec29bc0a264">theta</a> = theta        
<a name="l01588"></a>01588                       
<a name="l01589"></a>01589         <span class="comment"># write the output</span>
<a name="l01590"></a>01590         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a3abdd027570f1e8ee8b5a1c3937ebecd">nc_year</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a>] = (self.cur_year)
<a name="l01591"></a>01591         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#af00e43fdcd78dc1cf989e23c83576cf6">nc_doy</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a>] = (self.cur_doy)
<a name="l01592"></a>01592         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a326ef99eeb5d1b35707081c5d556ccb7">nc_sm</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>,:,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a>+1] = theta
<a name="l01593"></a>01593         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a45ec029f327b4811451d3db36b93c9a6">nc_recharge</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a>] = recharge_day
<a name="l01594"></a>01594         self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a919d83b479e4087a89d7216f8570b55b">nc_aet</a>[self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a8c8bc13a21b1982842d69007c295a305" title="run the GLUE ###########################">ens</a>,self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a>] = aet_day
<a name="l01595"></a>01595         
<a name="l01596"></a>01596         <span class="comment"># print progress</span>
<a name="l01597"></a>01597         <span class="keywordflow">if</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a> == int(0.25*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a64367da04fa6c5ab496748cf7bbf621c">max_t</a>):
<a name="l01598"></a>01598             output_message = <span class="stringliteral">&#39;25 % completed&#39;</span>
<a name="l01599"></a>01599             self._colored_output(output_message, 32)
<a name="l01600"></a>01600         
<a name="l01601"></a>01601         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a> == int(0.5*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a64367da04fa6c5ab496748cf7bbf621c">max_t</a>):
<a name="l01602"></a>01602             output_message = <span class="stringliteral">&#39;50 % completed&#39;</span>
<a name="l01603"></a>01603             self._colored_output(output_message, 32)
<a name="l01604"></a>01604         
<a name="l01605"></a>01605         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a> == int(0.75*self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a64367da04fa6c5ab496748cf7bbf621c">max_t</a>):
<a name="l01606"></a>01606             output_message = <span class="stringliteral">&#39;75 % completed&#39;</span>
<a name="l01607"></a>01607             self._colored_output(output_message, 32)
<a name="l01608"></a>01608         
<a name="l01609"></a>01609         <span class="keywordflow">elif</span> self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#aed5d7fb5129f294ac7a82f218dda4de0" title="run the model ########################">t</a> == self.<a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html#a64367da04fa6c5ab496748cf7bbf621c">max_t</a>-1:
<a name="l01610"></a>01610             output_message = <span class="stringliteral">&#39;100 % completed&#39;</span>
<a name="l01611"></a>01611             self._colored_output(output_message, 32)
<a name="l01612"></a>01612         <span class="comment">#print self.t</span>
<a name="l01613"></a>01613         
<a name="l01614"></a>01614 <span class="keywordflow">if</span> __name__==<span class="stringliteral">&#39;__main__&#39;</span>:
<a name="l01615"></a>01615      
<a name="l01616"></a>01616     <span class="comment">#maddur = RICHARDS_1D(&#39;/home/tomer/richards/input/maddur.xls&#39;)</span>
<a name="l01617"></a>01617     <span class="comment">#output_file = nc.NetCDFFile(maddur.ofile_name, &#39;r&#39;)</span>
<a name="l01618"></a>01618     <span class="comment">#print output_file.variables</span>
<a name="l01619"></a>01619     <span class="comment">#foo = output_file.variables[&#39;sm&#39;]</span>
<a name="l01620"></a>01620     <span class="comment">#theta= foo.getValue()</span>
<a name="l01621"></a>01621     <span class="comment">#print theta[:,-2]</span>
<a name="l01622"></a>01622     <span class="comment">#print theta[:,-1]</span>
<a name="l01623"></a>01623     <span class="comment">#plt.plot(theta[:,-1]); plt.plot(theta[:,-2]); plt.show()</span>
<a name="l01624"></a>01624     
<a name="l01625"></a>01625     <span class="comment">#maddur_ens = RICHARDS_1D_ENKF(&#39;/home/tomer/richards/input/maddur_ens.xls&#39;)</span>
<a name="l01626"></a>01626     <span class="comment">#output_file = nc.NetCDFFile(maddur_ens.ofile_name, &#39;r&#39;)</span>
<a name="l01627"></a>01627     <span class="comment">#foo = output_file.variables[&#39;sm&#39;][:,0,:]</span>
<a name="l01628"></a>01628     
<a name="l01629"></a><a class="code" href="namespaceambhas_1_1richards.html#acccd073afc15e3cb67349bc7d34092f3">01629</a>     maddur_glue = <a class="code" href="classambhas_1_1richards_1_1RICHARDS__1D__GLUE.html">RICHARDS_1D_GLUE</a>(<span class="stringliteral">&#39;/home/tomer/richards/input/maddur_glue_57.xls&#39;</span>)
<a name="l01630"></a>01630     
<a name="l01631"></a>01631     
</pre></div></div><!-- contents -->
</div>
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark">&#160;</span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark">&#160;</span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark">&#160;</span>Namespaces</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark">&#160;</span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark">&#160;</span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark">&#160;</span>Variables</a></div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

  <div id="nav-path" class="navpath">
    <ul>
      <li class="navelem"><a class="el" href="richards_8py.html">richards.py</a>      </li>

    <li class="footer">Generated on Sat Jul 21 2012 12:26:08 for AMBHAS by
    <a href="http://www.doxygen.org/index.html">
    <img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.6.1 </li>
   </ul>
 </div>


</body>
</html>
